#include "GridField.h"
#include "RenderMiscs.h"
#include "DecoLogger.h"
#include "ConfigManager.h"
#include "tetVisualizer.h"
#include "ParticleSys.h"
#include "Scene.h"

void GridField::CopyFrom(const GridField& rhs)
{
	DestroyAllArrays();
	scnBox = rhs.scnBox;

	numXElement = rhs.numXElement;
	numYElement = rhs.numYElement;
	numZElement = rhs.numZElement;
	stepSize = rhs.stepSize;

	numXVolume = rhs.numXVolume;
	numYVolume = rhs.numYVolume;
	numZVolume = rhs.numZVolume;
	stepSizeVol = rhs.stepSizeVol;

	if (numXElement)
	{
		SetNumElement(numXElement, numYElement, numZElement);
		CopyTo3DArray<DOUBLE>(u, rhs.u, (numXElement + 1), numYElement, numZElement);
		CopyTo3DArray<DOUBLE>(v, rhs.v, numXElement, (numYElement + 1), numZElement);
		CopyTo3DArray<DOUBLE>(w, rhs.w, numXElement, numYElement, (numZElement + 1));
		CopyTo3DArray<DOUBLE>(p, rhs.p, numXElement, numYElement, numZElement);
		CopyTo3DArray<DOUBLE>(phi, rhs.phi, numXElement, numYElement, numZElement);
		CopyTo3DArray<GridType>(type, rhs.type, numXElement, numYElement, numZElement);
//  [4/16/2008 HATEVOL] ADDED start
		CopyTo3DArray<GridType>(fastType, rhs.fastType, numXElement+2, numYElement+2, numZElement+2);
//  [4/16/2008 HATEVOL] ADDED end
		CopyTo3DArray<BOOL>(inDomain, rhs.inDomain, numXElement, numYElement, numZElement);
		CopyTo3DArray<BOOL>(prevInDomain, rhs.prevInDomain, numXElement, numYElement, numZElement);
		CopyTo3DArray<vector3>(pTheta, rhs.pTheta, numXElement, numYElement, numZElement);
		CopyTo3DArray<vector3>(mTheta, rhs.mTheta, numXElement, numYElement, numZElement);
		CopyTo3DArray<vector3>(mVelocity, rhs.mVelocity, numXElement, numYElement, numZElement);
		CopyTo3DArray<DOUBLE>(mVorticity, rhs.mVorticity, numXElement, numYElement, numZElement);
		CopyTo3DArray<DecoSolidObject*>(pSolid, rhs.pSolid, numXElement, numYElement, numZElement);
#ifdef _SURFACE_TENSION
		CopyTo3DArray<vector3>(pPj, rhs.pPj, numXElement, numYElement, numZElement);
		CopyTo3DArray<vector3>(mPj, rhs.mPj, numXElement, numYElement, numZElement);
#endif

	}
}

void GridField::CopyVelocityFrom(const GridField& rhs, UINT numX, UINT numY, UINT numZ)
{
	scnBox = rhs.scnBox;
	assert(numX <= rhs.numXElement);
	assert(numY <= rhs.numYElement);
	assert(numZ <= rhs.numZElement);
	assert(numX <= numXElement);
	assert(numY <= numYElement);
	assert(numZ <= numZElement);
	stepSize = rhs.stepSize;

	numXVolume = rhs.numXVolume;
	numYVolume = rhs.numYVolume;
	numZVolume = rhs.numZVolume;
	stepSizeVol = rhs.stepSizeVol;

	//DestroyAllArrays();
	if (numXElement)
	{
		//u = Malloc3DArray<DOUBLE>(numXElement + 1, numYElement, numZElement, TRUE);
		//v = Malloc3DArray<DOUBLE>(numXElement, numYElement + 1, numZElement, TRUE);
		//w = Malloc3DArray<DOUBLE>(numXElement, numYElement, numZElement + 1, TRUE);
		CopyTo3DArray<DOUBLE>(u, rhs.u, (numX + 1), numY, numZ);
		CopyTo3DArray<DOUBLE>(v, rhs.v, numX, (numY + 1), numZ);
		CopyTo3DArray<DOUBLE>(w, rhs.w, numX, numY, (numZ + 1));
		CopyTo3DArray<BOOL>(inDomain, rhs.inDomain, numX, numY, numZ);		
		CopyTo3DArray<GridType>(type, rhs.type, numXElement, numYElement, numZElement);
		CopyTo3DArray<GridType>(fastType, rhs.fastType, numXElement+2, numYElement+2, numZElement+2);
	}
}

void GridField::SetNumElement (UINT numX, UINT numY, UINT numZ, BOOL bFullAlloc)
{
	DestroyAllArrays();
	numXElement = numX;
	numYElement = numY;
	numZElement = numZ;

	stepSize = 2 * scnBox.GetExtent().x / numXElement;

	volxelPrec = 3;
	numXVolume = numXElement * volxelPrec;
	numYVolume = numYElement * volxelPrec;
	numZVolume = numZElement * volxelPrec;
	stepSizeVol = stepSize / volxelPrec;

	u = Malloc3DArray<DOUBLE>(numXElement + 1, numYElement, numZElement, TRUE);
	v = Malloc3DArray<DOUBLE>(numXElement, numYElement + 1, numZElement, TRUE);
	w = Malloc3DArray<DOUBLE>(numXElement, numYElement, numZElement + 1, TRUE);
	inDomain = Malloc3DArray<BOOL>(numXElement, numYElement, numZElement, TRUE);
	type = Malloc3DArray<GridType>(numXElement, numYElement, numZElement, TRUE);
	//  [4/16/2008 HATEVOL] ADDED start
	fastType = Malloc3DArray<GridType>(numXElement+2, numYElement+2, numZElement+2, TRUE);
	//  [4/16/2008 HATEVOL] ADDED end

	if (bFullAlloc)
	{
		p = Malloc3DArray<DOUBLE>(numXElement, numYElement, numZElement, TRUE);
		phi = Malloc3DArray<DOUBLE>(numXElement, numYElement, numZElement, TRUE);
		uFluidVol = Malloc3DArray<FLOAT>(numXElement + 1, numYElement, numZElement, TRUE);
		vFluidVol = Malloc3DArray<FLOAT>(numXElement, numYElement + 1, numZElement, TRUE);
		wFluidVol = Malloc3DArray<FLOAT>(numXElement, numYElement, numZElement + 1, TRUE);
		cFluidVol = Malloc3DArray<FLOAT>(numXElement, numYElement, numZElement, TRUE);
		tmpInDomain = Malloc3DArray<BOOL>(numXElement + 1, numYElement + 1, numZElement + 1);
		prevInDomain = Malloc3DArray<BOOL>(numXElement, numYElement, numZElement, TRUE);
		pTheta = Malloc3DArray<vector3>(numXElement, numYElement, numZElement, TRUE);
		mTheta = Malloc3DArray<vector3>(numXElement, numYElement, numZElement, TRUE);
		//testInDomain = Malloc3DArray<BOOL>(numXElement, numYElement, numZElement, TRUE);
		mVelocity = Malloc3DArray<vector3>(numXElement, numYElement, numZElement, TRUE);
		mVorticity = Malloc3DArray<DOUBLE>(numXElement, numYElement, numZElement, TRUE);
		pSolid = Malloc3DArray<DecoSolidObject*>(numXElement, numYElement, numZElement, TRUE);
	#ifdef _SURFACE_TENSION
		pPj = Malloc3DArray<vector3>(numXElement, numYElement, numZElement, TRUE);
		mPj = Malloc3DArray<vector3>(numXElement, numYElement, numZElement, TRUE);
	#endif
		for (UINT i = 0; i < numX; i++)
		{
			for (UINT j = 0; j < numY; j++)
			{
				for (UINT k = 0; k < numZ; k++)
				{
					phi[i][j][k] = MAX_DOUBLE;
				}
			}
		}
	}

}
void GridField::DestroyAllArrays()
{
	if (p)
		Free3DArray(p, numXElement, numYElement, numZElement);
	if (phi)
		Free3DArray(phi, numXElement, numYElement, numZElement);
	if (u)
		Free3DArray(u, numXElement + 1, numYElement, numZElement);
	if (v)
		Free3DArray(v, numXElement, numYElement + 1, numZElement);
	if (w)
		Free3DArray(w, numXElement, numYElement, numZElement + 1);
	if (uFluidVol)
		Free3DArray(uFluidVol, numXElement + 1, numYElement, numZElement);
	if (vFluidVol)
		Free3DArray(vFluidVol, numXElement, numYElement + 1, numZElement);
	if (wFluidVol)
		Free3DArray(wFluidVol, numXElement, numYElement, numZElement + 1);
	if (cFluidVol)
		Free3DArray(cFluidVol, numXElement, numYElement, numZElement);
	if (type)
		Free3DArray(type, numXElement, numYElement, numZElement);
//  [4/16/2008 HATEVOL] ADDED start
	if (fastType)
		Free3DArray(fastType, numXElement+2, numYElement+2, numZElement+2);
//  [4/16/2008 HATEVOL] ADDED end
	if (inDomain)
		Free3DArray(inDomain, numXElement, numYElement, numZElement);
	if (prevInDomain)
		Free3DArray(prevInDomain, numXElement, numYElement, numZElement);
	if (pTheta)
		Free3DArray(pTheta, numXElement, numYElement, numZElement);
	if (mTheta)
		Free3DArray(mTheta, numXElement, numYElement, numZElement);
	//if (testInDomain)
	//	Free3DArray(testInDomain, numXElement, numYElement, numZElement);
	if (tmpInDomain)
		Free3DArray(tmpInDomain, numXElement + 1, numYElement + 1, numZElement + 1);
	if (mVelocity)
		Free3DArray(mVelocity, numXElement, numYElement, numZElement);
	if (mVorticity)
		Free3DArray(mVorticity, numXElement, numYElement, numZElement);
	if (pSolid)
		Free3DArray(pSolid, numXElement, numYElement, numZElement);
#ifdef _SURFACE_TENSION
	if (pPj)
		Free3DArray(pPj, numXElement, numYElement, numZElement);
	if (mPj)
		Free3DArray(mPj, numXElement, numYElement, numZElement);
#endif
}

vector3 GridField::GetGridCenter(UINT xIndex, UINT yIndex, UINT zIndex) const
{
	vector3 scnBoxCenter;
	vector3 scnBoxExt;
	scnBox.GetCenterAndExtents(scnBoxCenter, scnBoxExt);
	vector3 scnBoxMin = scnBoxCenter - scnBoxExt;
	return vector3((xIndex + 0.5) * stepSize, (yIndex + 0.5) * stepSize, (zIndex + 0.5) * stepSize) + scnBoxMin;
}

vector3 GridField::GetSubGridCenter(UINT xIndex, UINT yIndex, UINT zIndex) const
{
	vector3 scnBoxCenter;
	vector3 scnBoxExt;
	scnBox.GetCenterAndExtents(scnBoxCenter, scnBoxExt);
	vector3 scnBoxMin = scnBoxCenter - scnBoxExt;
	return vector3((xIndex + 0.5) * stepSizeVol, (yIndex + 0.5) * stepSizeVol, (zIndex + 0.5) * stepSizeVol) + scnBoxMin;
}

Box GridField::GetGridBox(UINT xIndex, UINT yIndex, UINT zIndex) const
{
	vector3 center = GetGridCenter(xIndex, yIndex, zIndex);
	return Box(center - vector3(stepSize, stepSize, stepSize) / 2, center + vector3(stepSize, stepSize, stepSize) / 2);
}


void GridField::SetCellType(UINT xIndex, UINT yIndex, UINT zIndex, GridType newType)
{
	if (newType == GT_Air || newType == GT_Solid)
	{
		inDomain[xIndex][yIndex][zIndex] = FALSE;
		p[xIndex][yIndex][zIndex] = 0;
	}
	if (newType == GT_Surface || newType == GT_FLuid)
	{
		inDomain[xIndex][yIndex][zIndex] = TRUE;
	}

	type[xIndex][yIndex][zIndex] = newType;
	fastType[xIndex + 1][yIndex + 1][zIndex + 1] = newType;
}

void GridField::ClearTypeFlagTo(GridType newType)
{
//  [4/16/2008 HATEVOL] ADDED start
	for ( INT j = 0; j < numYElement+2; j++ ){
		for ( INT k = 0; k < numZElement+2; k++ ){
			fastType[0][j][k] = GT_NULL;
			fastType[numXElement+1][j][k] = GT_NULL;
		}
	}
	for (INT i = 0; i < numXElement; i++)
	{
		for ( INT k = 0; k < numZElement+2; k++ ){
			fastType[i+1][0][k] = GT_NULL;
			fastType[i+1][numYElement+1][k] = GT_NULL;
		}
		for (INT j = 0; j < numYElement; j++)
		{
			fastType[i+1][j+1][0] = GT_NULL;
			fastType[i+1][j+1][numZElement+1] = GT_NULL;
			for (INT k = 0; k < numZElement; k++)
			{
				if (newType == GT_Air)
				{
					inDomain[i][j][k] = FALSE;
					p[i][j][k] = 0;
				}
				if (newType == GT_Surface || newType == GT_FLuid)
				{
					inDomain[i][j][k] = TRUE;
				}
				type[i][j][k] = newType;
				fastType[i+1][j+1][k+1] = newType;
			}
		}
	}
//  [4/16/2008 HATEVOL] ADDED end
}

void GridField::Render(DecoRenderInterface* RI, DecoLight** Lights, UINT numEffectiveLights)
{
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	glRotatef(90.0f, 1, 0, 0);
	matrix44 mat;
	glGetDoublev(GL_MODELVIEW_MATRIX, (double*)&mat);
	glLoadIdentity();

	matrix44 mat2;
	RI->GetTransform(TT_LocalToWorld, mat2);
	RI->SetTransform(TT_LocalToWorld, mat * mat2);

	DecoRenderMisc::GetSingleton()->DrawBox(NULL, scnBox, DecoColor(255, 255, 255));

	UINT numElement = numXElement * numYElement * numZElement;
	if (!numZElement) return;
	for (INT xIndex = 0; xIndex < numXElement; xIndex++)
	{
		for (INT yIndex = 0; yIndex < numYElement; yIndex++)
		{
			for (INT zIndex = 0; zIndex < numZElement; zIndex++)
			{
//				if (zIndex != numZElement / 2 && zIndex != numZElement / 2 - 1)
//					continue;
				vector3 center = GetGridCenter(xIndex, yIndex, zIndex);
				Box box(center - vector3(stepSize, stepSize, stepSize) / 2, center + vector3(stepSize, stepSize, stepSize) / 2);
				DecoColor color;
				if (inDomain[xIndex][yIndex][zIndex])
				{
					color = DecoColor(255, 0, 255);
					//DecoRenderMisc::GetSingleton()->DrawBox(NULL, box, color);
				}
				switch (type[xIndex][yIndex][zIndex])
				{
				case GT_Air:
					color = DecoColor(255, 255, 255);
					break;
				case GT_FLuid:
					{
						color = DecoColor(255, 0, 255);
	//					DecoRenderMisc::GetSingleton()->DrawBox(NULL, box, color);
						break; 			
					}
				case GT_Solid:
					color = DecoColor(0, 255, 0);	
					DecoRenderMisc::GetSingleton()->DrawBox(NULL, box, color);
					break;
				case GT_Surface:
					color = DecoColor(255, 0, 255);
					//DecoRenderMisc::GetSingleton()->DrawBox(NULL, box, color);
					break;

				default:
					assert(0);
				}

				//DOUBLE tmpPhi = abs(phi[xIndex][yIndex][zIndex]);


				//BYTE cl = static_cast<BYTE>(Clamp<DOUBLE>(255 - tmpPhi * 15, 0, 254.99)) ;
				//color = DecoColor(cl, cl, cl);
				//vector3 gridCenter = GetGridCenter(xIndex, yIndex, zIndex);
				//DecoRenderMisc::GetSingleton()->DrawPointList(NULL, &gridCenter, 1, color, 2);

				//DOUBLE lineProp = 2 * stepSize;
				//DOUBLE velocityU = (u[xIndex][yIndex][zIndex] + u[xIndex + 1][yIndex][zIndex]) / 2;
				//DOUBLE velocityV = (v[xIndex][yIndex][zIndex] + v[xIndex][yIndex + 1][zIndex]) / 2;
				//DOUBLE velocityW = (w[xIndex][yIndex][zIndex] + w[xIndex][yIndex][zIndex + 1]) / 2;
				//vector3 velocity(velocityU, velocityV, velocityW);
				//if (velocity.lengthSqr() > EPSILON_FLOAT)
				//{
				//	vector3 lineStart = GetGridCenter(xIndex, yIndex, zIndex);
				//	DecoRenderMisc::GetSingleton()->DrawPointList(NULL, &lineStart, 1, DecoColor(0xffffffff), 2);
				//	vector3 lineEnd = lineStart + lineProp * velocity;
				//	DecoRenderMisc::GetSingleton()->DrawLine(lineStart, lineEnd, DecoColor(0xffffffff));
				//}

			}
		}

	}
	RI->SetTransform(TT_LocalToWorld, mat2);
}


void GridField::updateVoxelizeType(BYTE*** voxElements, INT numX, INT numY, INT numZ, GridType gridType, vector<Index3>* updateDomain, INT id)
{	
	if (updateDomain)
		updateDomain->clear();

	const BYTE mask = 1 << id;

	for (INT i = 0; i < numXElement; i++)
	{
		for (INT j = 0; j < numYElement; j++)
		{
			for (INT k = 0; k < numZElement; k++)
			{
				if (type[i][j][k] != GT_Solid)
				{
					type[i][j][k] = GT_Air;
				    fastType[i+1][j+1][k+1] = GT_Air;
				}

				BOOL newlyAdded = FALSE;

				if (type[i][j][k] != gridType)
				{
					BOOL hasSolid = FALSE, hasFluid = FALSE;
					INT subGridCount = volxelPrec * volxelPrec * volxelPrec;
					INT count = 0;
					for (INT x = 0; x < volxelPrec; x++)
					{
						for (INT y = 0; y < volxelPrec; y++)
						{
							for (INT z = 0; z < volxelPrec; z++)
							{
								if((voxElements[i * volxelPrec + x][j * volxelPrec + y][k * volxelPrec + z] & mask) != 0)
									count++;
							}
						}
					}

					if(count / (FLOAT)subGridCount >= 0.5f)
					{
						type[i][j][k] = gridType;
						fastType[i+1][j+1][k+1] = gridType;
						newlyAdded = TRUE;
					}
// 					else if(hasSolid && hasFluid)
// 					{
// 						newlyAdded = TRUE;
// 					}
				}

				if (newlyAdded && updateDomain)
					updateDomain->push_back(Index3(i, j, k));
			}
		}
	}

}

void GridField::VoxelizeSprayer(BYTE*** sprayerElements, INT numX, INT numY, INT numZ)
{
	m_compDomain.clear();
	updateVoxelizeType(sprayerElements, numX, numY, numZ, GT_Sprayer, &m_compDomain);
	INT domainSize = static_cast<INT>(m_compDomain.size());
	vector3 sprayerVelocity(0, 0, 0);
	DecoConfig::GetSingleton()->GetVector3("Sprayer", "velocity", sprayerVelocity);
	for (INT it = 0; it < domainSize; it++)
	{
		INT i = m_compDomain[it].m_i;
		INT j = m_compDomain[it].m_j;
		INT k = m_compDomain[it].m_k;	
		u[i][j][k] = u[i + 1][j][k] = sprayerVelocity.x;
		v[i][j][k] = v[i][j + 1][k] = sprayerVelocity.y;
		w[i][j][k] = w[i][j][k + 1] = sprayerVelocity.z;
	}
}
const vector<Index3>& GridField::updateSolidType(DecoSceneObject& solid, BYTE*** solidElements, INT numX, INT numY, INT numZ, INT id)
{
	m_compDomain.clear();
	m_validU.clear();
	m_validV.clear();
	m_validW.clear();

	updateVoxelizeType(solidElements, numX, numY, numZ, GT_Solid, &m_compDomain, id);
	//solid.SetAssociateGridIndex(m_compDomain);
	assert(solid.GetType() == SolidT || solid.GetType() == ArticulatedT);
	DecoSolidObject* solidObj = static_cast<DecoSolidObject*>(&solid);
	INT domainSize = static_cast<INT>(m_compDomain.size());
	for (INT it = 0; it < domainSize; ++it)
	{		
		INT i = m_compDomain[it].m_i;
		INT j = m_compDomain[it].m_j;
		INT k = m_compDomain[it].m_k;
		
		pSolid[i][j][k] = solidObj;
	}
	Zero3DArray(tmpInDomain, numXElement + 1, numYElement + 1, numZElement + 1);
	for (INT it = 0; it < domainSize; it++)
	{
		INT i = m_compDomain[it].m_i;
		INT j = m_compDomain[it].m_j;
		INT k = m_compDomain[it].m_k;
		if (!tmpInDomain[i][j][k])
		{
			m_validU.push_back(Index3(i, j, k));
			tmpInDomain[i][j][k] = TRUE;
		}
		if (!tmpInDomain[i + 1][j][k])
		{
			m_validU.push_back(Index3(i + 1, j, k));
			tmpInDomain[i + 1][j][k] = TRUE;
		}
	}
	Zero3DArray(tmpInDomain, numXElement + 1, numYElement + 1, numZElement + 1);
	for (INT it = 0; it < domainSize; it++)
	{
		INT i = m_compDomain[it].m_i;
		INT j = m_compDomain[it].m_j;
		INT k = m_compDomain[it].m_k;	
		if (!tmpInDomain[i][j][k])
		{
			m_validV.push_back(Index3(i, j, k));
			tmpInDomain[i][j][k] = TRUE;
		}
		if (!tmpInDomain[i + 1][j][k])
		{
			m_validV.push_back(Index3(i, j + 1, k));
			tmpInDomain[i][j + 1][k] = TRUE;
		}
	}

	Zero3DArray(tmpInDomain, numXElement + 1, numYElement + 1, numZElement + 1);
	for (INT it = 0; it < domainSize; it++)
	{
		INT i = m_compDomain[it].m_i;
		INT j = m_compDomain[it].m_j;
		INT k = m_compDomain[it].m_k;	
		if (!tmpInDomain[i][j][k])
		{
			m_validW.push_back(Index3(i, j, k));
			tmpInDomain[i][j][k] = TRUE;
		}
		if (!tmpInDomain[i + 1][j][k])
		{
			m_validW.push_back(Index3(i, j, k + 1));
			tmpInDomain[i][j][k + 1] = TRUE;
		}
	}
	return m_compDomain;
}

void GridField::updateFluidSolidRatio(BYTE*** solidElements, INT numX, INT numY, INT numZ)
{
	INT xRatio = numX / numXElement;
	INT yRatio = numY / numYElement;
	INT zRatio = numZ / numZElement;
	assert(!(xRatio % 2));
	assert(!(yRatio % 2));
	assert(!(zRatio % 2));
	if (!(xRatio == yRatio && yRatio == zRatio))
	{
		(*DecoLogger::GetSingleton()) << "Warning: updateFluidSolidRatio failed because xRatio yRatio and zRation differ.\n";
	}
	//for (INT i = 0; i < numXElement + 1; i++)
	//{
	//	for (INT j = 0; j < numYElement; j++)
	//	{
	//		for (INT k = 0; k < numZElement; k++)
	//		{
	//			if (!isSolidCell(i - 1, j, k) && !isSolidCell(i, j, k))
	//				continue;
	INT validUSize = static_cast<INT>(m_validU.size());
	for (INT it = 0; it < validUSize; it++)
	{
		INT i = m_validU[it].m_i;
		INT j = m_validU[it].m_j;
		INT k = m_validU[it].m_k;	
		if (xRatio > 0)
		{
			INT xEnd = ((i == numXElement) ? (numXElement * xRatio - 1) : (i * xRatio));
			INT xStart = ((i == 0) ? 0 : (xEnd - xRatio));
			//INT xStart = (i == 0 ? 0 : (i * xRatio - xRatio / 2));
			INT yStart = j * yRatio;
			INT zStart = k * zRatio;
			//INT xEnd = (i == numXElement ? xStart + xRatio / 2 : xStart + xRatio);
			INT yEnd = yStart + yRatio;
			INT zEnd = zStart + zRatio;
			INT solidCount = 0;
			INT totalCount = 0;
			for (INT iSub = xStart; iSub <= xEnd; iSub++)
			{
				for (INT jSub = yStart; jSub < yEnd; jSub++)
				{
					for (INT kSub = zStart; kSub < zEnd; kSub++)
					{
						totalCount++;
						if (solidElements[iSub][jSub][kSub] == 1)
							solidCount++;
					}
				}
			}
			if (totalCount)
				uFluidVol[i][j][k] = 1 - (FLOAT)solidCount / totalCount;
			else
				uFluidVol[i][j][k] = 1;
		}
		else
		{
			uFluidVol[i][j][k] = 0;
		}
		//	}
		//}
	}
	//for (INT i = 0; i < numXElement; i++)
	//{
	//	for (INT j = 0; j < numYElement + 1; j++)
	//	{
	//		for (INT k = 0; k < numZElement; k++)
	//		{
				//if (!isSolidCell(i, j - 1, k) && !isSolidCell(i, j, k))
				//	continue;
	INT validVSize = static_cast<INT>(m_validV.size());
	for (INT it = 0; it < validVSize; it++)
	{
		INT i = m_validV[it].m_i;
		INT j = m_validV[it].m_j;
		INT k = m_validV[it].m_k;	
		if (yRatio > 0)
		{
			INT yEnd = ((j == numYElement) ? (numYElement * yRatio - 1) : (j * yRatio));
			INT yStart = ((j == 0) ? 0 : (yEnd - yRatio));
			INT xStart = i * xRatio;
			//INT yStart = (j == 0 ? 0 : j * yRatio - yRatio / 2);
			INT zStart = k * zRatio;
			INT xEnd = xStart + xRatio;
			//INT yEnd = (j == numYElement ? yStart + yRatio / 2 : yStart + yRatio);
			INT zEnd = zStart + zRatio;
			INT solidCount = 0;
			INT totalCount = 0;
			for (INT iSub = xStart; iSub < xEnd; iSub++)
			{
				for (INT jSub = yStart; jSub <= yEnd; jSub++)
				{
					for (INT kSub = zStart; kSub < zEnd; kSub++)
					{
						totalCount++;
						if (solidElements[iSub][jSub][kSub] == 1)
							solidCount++;
					}
				}
			}
			if (totalCount)
				vFluidVol[i][j][k] = 1 - (FLOAT)solidCount / totalCount;
			else
				vFluidVol[i][j][k] = 1;
		}
		else
			vFluidVol[i][j][k] = 0;
		//	}
		//}
	}
	//for (INT i = 0; i < numXElement; i++)
	//{
	//	for (INT j = 0; j < numYElement; j++)
	//	{
	//		for (INT k = 0; k < numZElement + 1; k++)
	//		{
	//			if (!isSolidCell(i, j, k - 1) && !isSolidCell(i, j, k))
	//				continue;
	INT validWSize = static_cast<INT>(m_validW.size());
	for (INT it = 0; it < validWSize; it++)
	{
		INT i = m_validW[it].m_i;
		INT j = m_validW[it].m_j;
		INT k = m_validW[it].m_k;	
		if (zRatio > 0)
		{
			INT zEnd = ((k == numZElement) ? (numZElement * zRatio - 1) : (k * zRatio));
			INT zStart = ((k == 0) ? 0 : (zEnd - zRatio));
			INT xStart = i * xRatio;
			INT yStart = j * yRatio;
			//INT zStart = (k == 0 ? 0 : k * zRatio - zRatio / 2);
			INT xEnd = xStart + xRatio;
			INT yEnd = yStart + yRatio;
			//INT zEnd = (k == numZElement ? zStart + zRatio / 2 : zStart + zRatio);
			INT solidCount = 0;
			INT totalCount = 0;
			for (INT iSub = xStart; iSub < xEnd; iSub++)
			{
				for (INT jSub = yStart; jSub < yEnd; jSub++)
				{
					for (INT kSub = zStart; kSub <= zEnd; kSub++)
					{
						totalCount++;
						if (solidElements[iSub][jSub][kSub] == 1)
							solidCount++;
					}
				}
			}
			if (totalCount)
				wFluidVol[i][j][k] = 1 - (FLOAT)solidCount / totalCount;
			else
				wFluidVol[i][j][k] = 1;
		}
		else
			wFluidVol[i][j][k] = 0;
		//	}
		//}
	}
}

void GridField::updateSolidVelocity(DecoSceneObject& solid)
{
	DOUBLE totalDivergence = 0;
	INT uSize = static_cast<INT>(m_validU.size());
	for (INT it = 0; it < uSize; it++)
	{

		INT i = m_validU[it].m_i;
		INT j = m_validU[it].m_j;
		INT k = m_validU[it].m_k;

		vector3 pos = GetGridCenter(i, j, k) + vector3(-stepSize / 2, 0, 0);
		vector3 vel = solid.SampleVelocity(pos);
		u[i][j][k] = vel.x;
	}
	INT vSize = static_cast<INT>(m_validV.size());
	for (INT it = 0; it < vSize; it++)
	{
		INT i = m_validV[it].m_i;
		INT j = m_validV[it].m_j;
		INT k = m_validV[it].m_k;
		
		vector3 pos = GetGridCenter(i, j, k) + vector3(0, -stepSize / 2,  0);
		vector3 vel = solid.SampleVelocity(pos);
		v[i][j][k] = vel.y;
	}
	INT wSize = static_cast<INT>(m_validW.size());
	for (INT it = 0; it < wSize; it++)
	{
		INT i = m_validW[it].m_i;
		INT j = m_validW[it].m_j;
		INT k = m_validW[it].m_k;
		vector3 pos = GetGridCenter(i, j, k) + vector3(0, 0, -stepSize / 2);
		vector3 vel = solid.SampleVelocity(pos);
		w[i][j][k] = vel.z;
	}
	//checkSolidGridDivergence(solid);

}

void GridField::checkSolidGridDivergence(DecoSceneObject& solid)
{
	for (int i = 0; i < numXElement; ++i)
	{
		for (int j = 0; j < numYElement; ++j)
		{
			for (int k = 0; k < numZElement; ++k)
			{
				if (isSolidCell(i, j, k))
				{
					double divergence = u[i + 1][j][k] - u[i][j][k] + v[i][j + 1][k] - v[i][j][k] + w[i][j][k + 1] - w[i][j][k];
					if (abs(divergence) > 1e-2)
					{
						vector3 pos = GetGridCenter(i, j, k) + vector3(0, -stepSize / 2,  0);
						vector3 vel = solid.SampleVelocity(pos);

						pos = GetGridCenter(i, j, k) + vector3(-stepSize / 2, 0,  0);
						vel = solid.SampleVelocity(pos);


						(*DecoLogger::GetSingleton()) << "divergence of (" << i << "," << j << "," << k << ") is " << divergence << "\n";
					}
				}
			}
		}
	}
}

BOOL GridField::isSolidBoundaryCell(INT xIndex, INT yIndex, INT zIndex)
{
	if (xIndex < 0 || yIndex < 0 || zIndex < 0 || xIndex >= static_cast<INT>(numXElement) || yIndex >= static_cast<INT>(numYElement) || zIndex >= static_cast<INT>(numZElement))
		return FALSE;
	if (type[xIndex][yIndex][zIndex] == GT_Solid &&
	    (type[xIndex + 1][yIndex][zIndex] != GT_Solid || type[xIndex - 1][yIndex][zIndex] != GT_Solid ||
		 type[xIndex][yIndex + 1][zIndex] != GT_Solid || type[xIndex][yIndex - 1][zIndex] != GT_Solid ||
		 type[xIndex][yIndex][zIndex + 1] != GT_Solid || type[xIndex][yIndex][zIndex - 1] != GT_Solid ))
		 return TRUE;
	else
		return FALSE;

}
const vector<Index3>& GridField::VoxelizeSolid(DecoSceneObject& solid, BYTE*** solidElements, INT numX, INT numY, INT numZ)
{

	const vector<Index3>& assoIndices = updateSolidType(solid, solidElements, numX, numY, numZ, 0);

	updateSolidVelocity(solid);
	return assoIndices;
}

#ifdef _SLC

BOOL GridField::IsThinSheetCell(INT xIndex, INT yIndex, INT zIndex)
{
	return (!isFluidCell(xIndex, yIndex, zIndex) && 
		!isFluidCell(xIndex - 1, yIndex, zIndex) && !isFluidCell(xIndex + 1, yIndex, zIndex) &&
		!isFluidCell(xIndex, yIndex - 1, zIndex) && !isFluidCell(xIndex, yIndex + 1, zIndex) &&
		!isFluidCell(xIndex, yIndex, zIndex - 1) && !isFluidCell(xIndex, yIndex, zIndex + 1));
}
BOOL GridField::IsInvalidAirCell(INT i, INT j, INT k)
{
	if (!isAirCell(i, j, k))
		return FALSE;
	if (isSolidCell(i, j - 1, k) && isSolidCell(i - 1, j, k))
	{
		if (isFluidCell(i, j + 1, k) && isFluidCell(i + 1, j, k))
			return TRUE;
	}
	if (isSolidCell(i, j - 1, k) && isSolidCell(i + 1, j, k))
	{
		if (isFluidCell(i, j + 1, k) && isFluidCell(i - 1, j, k))
			return TRUE;
	}
	if (isSolidCell(i, j - 1, k) && isSolidCell(i, j, k - 1))
	{
		if (isFluidCell(i, j + 1, k) && isFluidCell(i, j, k + 1))
			return TRUE;
	}
	if (isSolidCell(i, j - 1, k) && isSolidCell(i, j, k + 1))
	{
		if (isFluidCell(i, j + 1, k) && isFluidCell(i, j, k - 1))
			return TRUE;
	}
	if (isSolidCell(i - 1, j, k) && isSolidCell(i, j, k - 1))
	{
		if (isFluidCell(i + 1, j, k) && isFluidCell(i, j, k + 1))
			return TRUE;
	}
	if (isSolidCell(i - 1, j, k) && isSolidCell(i, j, k + 1))
	{
		if (isFluidCell(i + 1, j, k) && isFluidCell(i, j, k - 1))
			return TRUE;
	}
	if (isSolidCell(i + 1, j, k) && isSolidCell(i, j, k - 1))
	{
		if (isFluidCell(i - 1, j, k) && isFluidCell(i, j, k + 1))
			return TRUE;
	}
	if (isSolidCell(i + 1, j, k) && isSolidCell(i, j, k + 1))
	{
		if (isFluidCell(i - 1, j, k) && isFluidCell(i, j, k - 1))
			return TRUE;
	}
	return FALSE;
}
#endif


BOOL GridField::isFluidCell(INT xIndex, INT yIndex, INT zIndex) const
{

	return ( fastType[xIndex+1][yIndex+1][zIndex+1] == GT_FLuid || fastType[xIndex+1][yIndex+1][zIndex+1] == GT_Surface );
}
BOOL GridField::isAirCell(INT xIndex, INT yIndex, INT zIndex) const
{

	return (fastType[xIndex+1][yIndex+1][zIndex+1] == GT_Air);
}

BOOL GridField::isSolidCell (INT xIndex, INT yIndex, INT zIndex) const
{
	return ( fastType[xIndex+1][yIndex+1][zIndex+1] == GT_Solid || fastType[xIndex+1][yIndex+1][zIndex+1] == GT_NULL );
}

BOOL GridField::isSurfaceCell(INT xIndex, INT yIndex, INT zIndex) const
{
	return fastType[xIndex+1][yIndex+1][zIndex+1] == GT_Surface;
}

BOOL GridField::isSprayerCell(INT xIndex, INT yIndex, INT zIndex) const
{
	return fastType[xIndex+1][yIndex+1][zIndex+1] == GT_Sprayer;
}


BOOL GridField::isBoundaryCell(INT xIndex, INT yIndex, INT zIndex) const
{
	if (isSolidCell(xIndex - 1, yIndex, zIndex) || isSolidCell(xIndex + 1, yIndex, zIndex) || 
		isSolidCell(xIndex, yIndex - 1, zIndex) || isSolidCell(xIndex, yIndex + 1, zIndex) ||
		isSolidCell(xIndex, yIndex, zIndex - 1) || isSolidCell(xIndex, yIndex, zIndex + 1)) 
		return TRUE;
	else if (!inDomain[xIndex - 1][yIndex][zIndex] || !inDomain[xIndex + 1][yIndex][zIndex] ||
		!inDomain[xIndex][yIndex - 1][zIndex] || !inDomain[xIndex][yIndex + 1][zIndex] ||
		!inDomain[xIndex][yIndex][zIndex - 1] || !inDomain[xIndex][yIndex][zIndex + 1])
		return TRUE;
	else if (type[xIndex][yIndex][zIndex] == GT_Surface)
		return TRUE;
	if (isSprayerCell(xIndex - 1, yIndex, zIndex) || isSprayerCell(xIndex + 1, yIndex, zIndex) || 
		isSprayerCell(xIndex, yIndex - 1, zIndex) || isSprayerCell(xIndex, yIndex + 1, zIndex) ||
		isSprayerCell(xIndex, yIndex, zIndex - 1) || isSprayerCell(xIndex, yIndex, zIndex + 1)) 
		return TRUE;
	else
		return FALSE;
}
BOOL GridField::isNearSolid(INT xIndex, INT yIndex, INT zIndex) const
{
	if (isSolidCell(xIndex - 1, yIndex, zIndex) || isSolidCell(xIndex + 1, yIndex, zIndex) || 
		isSolidCell(xIndex, yIndex - 1, zIndex) || isSolidCell(xIndex, yIndex + 1, zIndex) ||
		isSolidCell(xIndex, yIndex, zIndex - 1) || isSolidCell(xIndex, yIndex, zIndex + 1)) 
		return true;
	else
		return false;
}
BOOL GridField::isInDomain(INT xIndex, INT yIndex, INT zIndex) const
{
	if (xIndex < 0 || xIndex >= numXElement || yIndex < 0 || yIndex >= numYElement || zIndex < 0 || zIndex >= numZElement)
		return FALSE;
	return inDomain[xIndex][yIndex][zIndex];
}

vector3 GridField::fromLocalGridToWorld(const vector3& pos, INT xIndex, INT yIndex, INT zIndex) const
{
	assert (pos.x <= stepSize && pos.y <= stepSize && pos.z <= stepSize);
	vector3 result;
	vector3 Min = scnBox.GetMin();
	result.x = xIndex * stepSize + pos.x + Min.x;
	result.y = yIndex * stepSize + pos.y + Min.y;
	result.z = zIndex * stepSize + pos.z + Min.z;
	return result;
}
vector3 GridField::fromWorldToLocalGrid(const vector3& pos, UINT& xIndex, UINT& yIndex, UINT& zIndex) const
{
	const vector3& Min = scnBox.GetMin();
	DOUBLE invertSteSize = 1.0 / stepSize;
	DOUBLE smallEpsilon = 0.01 * EPSILON_FLOAT;
	vector3 offset(pos - Min);
	//DOUBLE tmpXIndex = (offset.x * invertSteSize);
	//DOUBLE tmpYIndex = (offset.y * invertSteSize);
	//DOUBLE tmpZIndex = (offset.z * invertSteSize);

	INT tmpXIndex = static_cast<INT>(offset.x * invertSteSize);
	INT tmpYIndex = static_cast<INT>(offset.y * invertSteSize);
	INT tmpZIndex = static_cast<INT>(offset.z * invertSteSize);

	vector3 result;

	if ( tmpXIndex < 0 ){
		xIndex = 0;
		result.x = 0;
	}
	else if ( tmpXIndex > numXElement - 1 ) {
		xIndex = numXElement - 1;
		result.x = stepSize - smallEpsilon;
	}
	else {
		xIndex = static_cast<INT>(tmpXIndex);
		result.x = offset.x - xIndex * stepSize;
	}
	if ( tmpYIndex < 0 ){
		yIndex = 0;
		result.y = 0;
	}
	else if ( tmpYIndex > numYElement - 1 ) {
		yIndex = numYElement - 1;
		result.y = stepSize - smallEpsilon;
	}
	else {
		yIndex = static_cast<INT>(tmpYIndex);
		result.y = offset.y - yIndex * stepSize;
	}
	if ( tmpZIndex < 0 ){
		zIndex = 0;
		result.z = 0;
	}
	else if ( tmpZIndex > numZElement - 1 ) {
		zIndex = numZElement - 1;
		result.z = stepSize - smallEpsilon;
	}
	else {
		zIndex = static_cast<INT>(tmpZIndex);
		result.z = offset.z - zIndex * stepSize;
	}

	return result;
}

DOUBLE GridField::SamplePhi(const vector3& pos) const
{
	return SampleScalerValue(pos, this->phi);
}

DOUBLE GridField::SamplePressure(const vector3& pos) const
{
	return SampleScalerValue(pos, this->p);
}

DOUBLE GridField::SampleScalerValue(const vector3& pos, DOUBLE*** scaler) const
{

	INT i, j, k;
	vector3 localPos = fromWorldToLocalGrid(pos, (UINT&)i, (UINT&)j, (UINT&)k);

	if (NearlyEquals(localPos, vector3(stepSize / 2, stepSize / 2, stepSize / 2), EPSILON_FLOAT))
		return (scaler[i][j][k]);

	DOUBLE halfStep = stepSize / 2;
	DOUBLE xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, xAlpha, yAlpha, zAlpha, result; 


	INT numX = numXElement;
	INT numY = numYElement;
	INT numZ = numZElement;

	if (localPos.x >= halfStep && localPos.y >= halfStep && localPos.z >= halfStep)
	{
		xAlpha = (localPos.x - halfStep) / stepSize;
		zAlpha = (localPos.z - halfStep) / stepSize;
		yAlpha = (localPos.y - halfStep) / stepSize;

		xValue1 = scaler[i][j][k];
		xValue2 = scaler[!isInDomain(i + 1, j, k) ? i : i + 1][j][k];
		xValue3 = scaler[i][j][(!isInDomain(i, j, k + 1)) ? k : k + 1];
		xValue4 = scaler[!isInDomain(i + 1, j, k) ? i : i + 1][j][(!isInDomain(i, j, k + 1)) ? k : k + 1];
		yValue1 = scaler[i][(!isInDomain(i, j + 1, k)) ? j : j + 1][k];
		yValue2 = scaler[!isInDomain(i + 1, j, k) ? i : i + 1][(!isInDomain(i, j + 1, k)) ? j : j + 1][k];
		yValue3 = scaler[i][(!isInDomain(i, j + 1, k)) ? j : j + 1][(!isInDomain(i, j, k + 1)) ? k : k + 1];
		yValue4 = scaler[!isInDomain(i + 1, j, k) ? i : i + 1][(!isInDomain(i, j + 1, k)) ? j : j + 1][(!isInDomain(i, j, k + 1)) ? k : k + 1];
	}
	else if (localPos.x >= halfStep && localPos.y >= halfStep && localPos.z < halfStep)
	{
		xAlpha = (localPos.x - halfStep) / stepSize;
		zAlpha = (halfStep - localPos.z) / stepSize;
		yAlpha = (localPos.y - halfStep) / stepSize;

		xValue1 = scaler[i][j][k];
		xValue2 = scaler[!isInDomain(i + 1, j, k) ? i : i + 1][j][k];
		xValue3 = scaler[i][j][(!isInDomain(i, j, k - 1)) ? k : k - 1];
		xValue4 = scaler[!isInDomain(i + 1, j, k) ? i : i + 1][j][(!isInDomain(i, j, k - 1)) ? k : k - 1];
		yValue1 = scaler[i][(!isInDomain(i, j + 1, k)) ? j : j + 1][k];
		yValue2 = scaler[!isInDomain(i + 1, j, k) ? i : i + 1][(!isInDomain(i, j + 1, k)) ? j : j + 1][k];
		yValue3 = scaler[i][(!isInDomain(i, j + 1, k)) ? j : j + 1][(!isInDomain(i, j, k - 1)) ? k : k - 1];
		yValue4 = scaler[!isInDomain(i + 1, j, k) ? i : i + 1][(!isInDomain(i, j + 1, k)) ? j : j + 1][(!isInDomain(i, j, k - 1)) ? k : k - 1];
	}
	else if (localPos.x < halfStep && localPos.y >= halfStep && localPos.z < halfStep)
	{
		xAlpha = (halfStep - localPos.x) / stepSize;
		zAlpha = (halfStep - localPos.z) / stepSize;
		yAlpha = (localPos.y - halfStep) / stepSize;

		xValue1 = scaler[i][j][k];
		xValue2 = scaler[!isInDomain(i - 1, j, k) ? i : i - 1][j][k];
		xValue3 = scaler[i][j][(!isInDomain(i, j, k - 1)) ? k : k - 1];
		xValue4 = scaler[!isInDomain(i - 1, j, k) ? i : i - 1][j][(!isInDomain(i, j, k - 1)) ? k : k - 1];
		yValue1 = scaler[i][(!isInDomain(i, j + 1, k)) ? j : j + 1][k];
		yValue2 = scaler[!isInDomain(i - 1, j, k) ? i : i - 1][(!isInDomain(i, j + 1, k)) ? j : j + 1][k];
		yValue3 = scaler[i][(!isInDomain(i, j + 1, k)) ? j : j + 1][(!isInDomain(i, j, k - 1)) ? k : k - 1];
		yValue4 = scaler[!isInDomain(i - 1, j, k) ? i : i - 1][(!isInDomain(i, j + 1, k)) ? j : j + 1][(!isInDomain(i, j, k - 1)) ? k : k - 1];
	}

	else if (localPos.x < halfStep && localPos.y >= halfStep && localPos.z >= halfStep)
	{
		xAlpha = (halfStep - localPos.x) / stepSize;
		zAlpha = (localPos.z - halfStep) / stepSize;
		yAlpha = (localPos.y - halfStep) / stepSize;

		xValue1 = scaler[i][j][k];
		xValue2 = scaler[!isInDomain(i - 1, j, k) ? i : i - 1][j][k];
		xValue3 = scaler[i][j][(!isInDomain(i, j, k + 1)) ? k : k + 1];
		xValue4 = scaler[!isInDomain(i - 1, j, k) ? i : i - 1][j][(!isInDomain(i, j, k + 1)) ? k : k + 1];
		yValue1 = scaler[i][(!isInDomain(i, j + 1, k)) ? j : j + 1][k];
		yValue2 = scaler[!isInDomain(i - 1, j, k) ? i : i - 1][(!isInDomain(i, j + 1, k)) ? j : j + 1][k];
		yValue3 = scaler[i][(!isInDomain(i, j + 1, k)) ? j : j + 1][(!isInDomain(i, j, k + 1)) ? k : k + 1];
		yValue4 = scaler[!isInDomain(i - 1, j, k) ? i : i - 1][(!isInDomain(i, j + 1, k)) ? j : j + 1][(!isInDomain(i, j, k + 1)) ? k : k + 1];
	}
	else if (localPos.x >= halfStep && localPos.y < halfStep && localPos.z >= halfStep)
	{
		xAlpha = (localPos.x - halfStep) / stepSize;
		zAlpha = (localPos.z - halfStep) / stepSize;
		yAlpha = (halfStep - localPos.y) / stepSize;

		xValue1 = scaler[i][j][k];
		xValue2 = scaler[!isInDomain(i + 1, j, k) ? i : i + 1][j][k];
		xValue3 = scaler[i][j][(!isInDomain(i, j, k + 1)) ? k : k + 1];
		xValue4 = scaler[!isInDomain(i + 1, j, k) ? i : i + 1][j][(!isInDomain(i, j, k + 1)) ? k : k + 1];
		yValue1 = scaler[i][(!isInDomain(i, j - 1, k)) ? j : j - 1][k];
		yValue2 = scaler[!isInDomain(i + 1, j, k) ? i : i + 1][(!isInDomain(i, j - 1, k)) ? j : j - 1][k];
		yValue3 = scaler[i][(!isInDomain(i, j - 1, k)) ? j : j - 1][(!isInDomain(i, j, k + 1)) ? k : k + 1];
		yValue4 = scaler[!isInDomain(i + 1, j, k) ? i : i + 1][(!isInDomain(i, j - 1, k)) ? j : j - 1][(!isInDomain(i, j, k + 1)) ? k : k + 1];
	}
	else if (localPos.x >= halfStep && localPos.y < halfStep && localPos.z < halfStep)
	{
		xAlpha = (localPos.x - halfStep) / stepSize;
		zAlpha = (halfStep - localPos.z) / stepSize;
		yAlpha = (halfStep - localPos.y) / stepSize;

		xValue1 = scaler[i][j][k];
		xValue2 = scaler[!isInDomain(i + 1, j, k) ? i : i + 1][j][k];
		xValue3 = scaler[i][j][(!isInDomain(i, j, k - 1)) ? k : k - 1];
		xValue4 = scaler[!isInDomain(i + 1, j, k) ? i : i + 1][j][(!isInDomain(i, j, k - 1)) ? k : k - 1];
		yValue1 = scaler[i][(!isInDomain(i, j - 1, k)) ? j : j - 1][k];
		yValue2 = scaler[!isInDomain(i + 1, j, k) ? i : i + 1][(!isInDomain(i, j - 1, k)) ? j : j - 1][k];
		yValue3 = scaler[i][(!isInDomain(i, j - 1, k)) ? j : j - 1][(!isInDomain(i, j, k - 1)) ? k : k - 1];
		yValue4 = scaler[!isInDomain(i + 1, j, k) ? i : i + 1][(!isInDomain(i, j - 1, k)) ? j : j - 1][(!isInDomain(i, j, k - 1)) ? k : k - 1];
	}
	else if (localPos.x < halfStep && localPos.y < halfStep && localPos.z < halfStep)
	{
		xAlpha = (halfStep - localPos.x) / stepSize;
		zAlpha = (halfStep - localPos.z) / stepSize;
		yAlpha = (halfStep - localPos.y) / stepSize;

		xValue1 = scaler[i][j][k];
		xValue2 = scaler[!isInDomain(i - 1, j, k) ? i : i - 1][j][k];
		xValue3 = scaler[i][j][(!isInDomain(i, j, k - 1)) ? k : k - 1];
		xValue4 = scaler[!isInDomain(i - 1, j, k) ? i : i - 1][j][(!isInDomain(i, j, k - 1)) ? k : k - 1];
		yValue1 = scaler[i][(!isInDomain(i, j - 1, k)) ? j : j - 1][k];
		yValue2 = scaler[!isInDomain(i - 1, j, k) ? i : i - 1][(!isInDomain(i, j - 1, k)) ? j : j - 1][k];
		yValue3 = scaler[i][(!isInDomain(i, j - 1, k)) ? j : j - 1][(!isInDomain(i, j, k - 1)) ? k : k - 1];
		yValue4 = scaler[!isInDomain(i - 1, j, k) ? i : i - 1][(!isInDomain(i, j - 1, k)) ? j : j - 1][(!isInDomain(i, j, k - 1)) ? k : k - 1];
	}
	else if (localPos.x < halfStep && localPos.y < halfStep && localPos.z >= halfStep)
	{
		xAlpha = (halfStep - localPos.x) / stepSize;
		zAlpha = (localPos.z - halfStep) / stepSize;
		yAlpha = (halfStep - localPos.y) / stepSize;

		xValue1 = scaler[i][j][k];
		xValue2 = scaler[!isInDomain(i - 1, j, k) ? i : i - 1][j][k];
		xValue3 = scaler[i][j][(!isInDomain(i, j, k + 1)) ? k : k + 1];
		xValue4 = scaler[!isInDomain(i - 1, j, k) ? i : i - 1][j][(!isInDomain(i, j, k + 1)) ? k : k + 1];
		yValue1 = scaler[i][(!isInDomain(i, j - 1, k)) ? j : j - 1][k];
		yValue2 = scaler[!isInDomain(i - 1, j, k) ? i : i - 1][(!isInDomain(i, j - 1, k)) ? j : j - 1][k];
		yValue3 = scaler[i][(!isInDomain(i, j - 1, k)) ? j : j - 1][(!isInDomain(i, j, k + 1)) ? k : k + 1];
		yValue4 = scaler[!isInDomain(i - 1, j, k) ? i : i - 1][(!isInDomain(i, j - 1, k)) ? j : j - 1][(!isInDomain(i, j, k + 1)) ? k : k + 1];
	}
	else
		assert(0);

	//assert(xAlpha >= 0 && xAlpha <= 1 && yAlpha >= 0 && yAlpha <= 1 && zAlpha >= 0 && zAlpha <= 1);
	result = TrilinearInterpolate<DOUBLE>(xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, xAlpha, zAlpha, yAlpha);
	return result;
}

//  [4/24/2008 HATEVOL] ADDED start

DOUBLE GridField::SampleUWithLocalGrid(const vector3& localPos, INT i, INT j, INT k, BOOL bPrintf /* = FALSE */) const
{
	DOUBLE resultU;
	const vector3 minPt = scnBox.GetMin();
	const vector3 maxPt = scnBox.GetMax();
	
	INT numX = numXElement;
	INT numY = numYElement;
	INT numZ = numZElement;
	//INT i = 0, j = 0, k = 0;

	//vector3 localPos = fromWorldToLocalGrid(pos, (UINT&)i, (UINT&)j, (UINT&)k);

	DOUBLE xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, xAlpha, yAlpha, zAlpha; 

	{
		DOUBLE invertStepSize = 1 / stepSize;
		xValue1 = u[i][j][k];
		xValue2 = u[i + 1][j][k];
		xAlpha = localPos.x * invertStepSize;
		zAlpha = signedAbs(0.5 - localPos.z * invertStepSize );
		yAlpha = signedAbs(0.5 - localPos.y * invertStepSize);

		DOUBLE halfStepSize = 0.5 * stepSize;
		
		BOOL isSolidIJKminus, isSolidIJKplus, isSolidIJminusK, isSolidIJplusK;

		if ( localPos.y >= halfStepSize )
		{
			if ( localPos.z >= halfStepSize )
			{
				isSolidIJKplus  = isSolidCell(i, j, k+1);
				isSolidIJplusK  = isSolidCell(i, j+1, k);
				xValue3 = u[i][j][isSolidIJKplus ? k : k + 1];
				yValue1 = u[i][isSolidIJplusK ? j : j + 1][k];
				yValue3 = u[i][isSolidIJplusK ? j : j + 1][isSolidIJKplus ? k : k + 1];
				xValue4 = u[i + 1][j][isSolidIJKplus ? k : k + 1];
				yValue2 = u[i + 1][isSolidIJplusK ? j : j + 1][k];
				yValue4 = u[i + 1][isSolidIJplusK ? j : j + 1][isSolidIJKplus ? k : k + 1];

			}
			else {
				isSolidIJKminus = isSolidCell(i, j, k-1);
				isSolidIJplusK  = isSolidCell(i, j+1, k);
				xValue3 = u[i][j][isSolidIJKminus ? k : k - 1];
				yValue1 = u[i][isSolidIJplusK ? j : j + 1][k];
				yValue3 = u[i][isSolidIJplusK ? j : j + 1][isSolidIJKminus ? k : k - 1];
				xValue4 = u[i + 1][j][isSolidIJKminus ? k : k - 1];
				yValue2 = u[i + 1][isSolidIJplusK ? j : j + 1][k];
				yValue4 = u[i + 1][isSolidIJplusK ? j : j + 1][isSolidIJKminus ? k : k - 1];
			}
		}
		else 
		{
			if ( localPos.z >= halfStepSize)
			{
				isSolidIJKplus  = isSolidCell(i, j, k+1);
				isSolidIJminusK = isSolidCell(i, j-1, k);
				xValue3 = u[i][j][isSolidIJKplus ? k : k + 1];
				yValue1 = u[i][isSolidIJminusK ? j : j - 1][k];
				yValue3 = u[i][isSolidIJminusK ? j : j - 1][isSolidIJKplus ? k : k + 1];
				xValue4 = u[i + 1][j][isSolidIJKplus ? k : k + 1]; 
				yValue2 = u[i + 1][isSolidIJminusK ? j : j - 1][k];
				yValue4 = u[i + 1][isSolidIJminusK ? j : j - 1][isSolidIJKplus ? k : k + 1];
			}
			else 
			{
				isSolidIJKminus = isSolidCell(i, j, k-1);
				isSolidIJminusK = isSolidCell(i, j-1, k);
				xValue3 = u[i][j][isSolidIJKminus ? k : k - 1];
				yValue1 = u[i][isSolidIJminusK ? j : j - 1][k];
				yValue3 = u[i][isSolidIJminusK ? j : j - 1][isSolidIJKminus ? k : k - 1];
				xValue4 = u[i + 1][j][isSolidIJKminus ? k : k - 1];
				yValue2 = u[i + 1][isSolidIJminusK ? j : j - 1][k];
				yValue4 = u[i + 1][isSolidIJminusK ? j : j - 1][isSolidIJKminus ? k : k - 1];

			}

		}
		resultU = TrilinearInterpolate<DOUBLE>(xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, xAlpha, zAlpha, yAlpha);
	}
	if (bPrintf)
	{
		(*DecoLogger::GetSingleton()) << "xValue1 :" << xValue1 << "\n";
		(*DecoLogger::GetSingleton()) << "xValue2 :" << xValue2 << "\n";
		(*DecoLogger::GetSingleton()) << "xValue3 :" << xValue3 << "\n";
		(*DecoLogger::GetSingleton()) << "xValue4 :" << xValue4 << "\n";
		(*DecoLogger::GetSingleton()) << "yValue1 :" << yValue1 << "\n";
		(*DecoLogger::GetSingleton()) << "yValue2 :" << yValue2 << "\n";
		(*DecoLogger::GetSingleton()) << "yValue3 :" << yValue3 << "\n";
		(*DecoLogger::GetSingleton()) << "yValue4 :" << yValue4 << "\n";
		(*DecoLogger::GetSingleton()) << "xAlpha :" << xAlpha << "\n";
		(*DecoLogger::GetSingleton()) << "yAlpha :" << yAlpha << "\n";
		(*DecoLogger::GetSingleton()) << "zAlpha :" << zAlpha << "\n";
		(*DecoLogger::GetSingleton()) << "resultU :" << resultU << "\n";
	}
	return resultU;
}

DOUBLE GridField::SampleVWithLocalGrid(const vector3& localPos, INT i, INT j, INT k) const
{
	DOUBLE resultV;
	const vector3 minPt = scnBox.GetMin();
	const vector3 maxPt = scnBox.GetMax();
	INT numX = numXElement;
	INT numY = numYElement;
	INT numZ = numZElement;
	//INT i = 0, j = 0, k = 0;

	//vector3 localPos = fromWorldToLocalGrid(pos, (UINT&)i, (UINT&)j, (UINT&)k);

	DOUBLE xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, xAlpha, yAlpha, zAlpha; 

	{
		DOUBLE invertStepSize = 1.0 / stepSize;
		xValue1 = v[i][j][k];
		xValue2 = v[i][j + 1][k];
		yAlpha = localPos.y * invertStepSize;
		zAlpha = signedAbs(0.5 - localPos.z * invertStepSize);
		xAlpha = signedAbs(0.5 - localPos.x * invertStepSize);
		
		DOUBLE halfStepSize = 0.5 * stepSize;
		BOOL isSolidIJKminus, isSolidIJKplus, isSolidIminusJK, isSolidIplusJK;
		if (localPos.x >= halfStepSize)
		{
			if ( localPos.z >= halfStepSize)
			{
				isSolidIJKplus = isSolidCell(i, j, k+1);
				isSolidIplusJK = isSolidCell(i+1, j, k);
				xValue3 = v[i][j][isSolidIJKplus ? k : k + 1];
				xValue4 = v[i][j + 1][isSolidIJKplus ? k : k + 1];
				yValue1 = v[isSolidIplusJK ? i : i + 1][j][k];
				yValue2 = v[isSolidIplusJK ? i : i + 1][j + 1][k];
				yValue3 = v[isSolidIplusJK ? i : i + 1][j][isSolidIJKplus ? k : k + 1];
				yValue4 = v[isSolidIplusJK ? i : i + 1][j + 1][isSolidIJKplus ? k : k + 1];
			}
			else
			{
				isSolidIJKminus = isSolidCell(i, j, k-1);
				isSolidIplusJK  = isSolidCell(i+1, j, k);
				xValue3 = v[i][j][isSolidIJKminus ? k : k - 1];
				xValue4 = v[i][j + 1][isSolidIJKminus ? k : k - 1];
				yValue1 = v[isSolidIplusJK ? i : i + 1][j][k];
				yValue2 = v[isSolidIplusJK ? i : i + 1][j + 1][k];
				yValue3 = v[isSolidIplusJK ? i : i + 1][j][isSolidIJKminus ? k : k - 1];
				yValue4 = v[isSolidIplusJK ? i : i + 1][j + 1][isSolidIJKminus ? k : k - 1];
			}
		}
		else 
		{
			if (localPos.z >= halfStepSize)
			{
				isSolidIJKplus  = isSolidCell(i, j, k+1);
				isSolidIminusJK = isSolidCell(i-1, j, k);
				xValue3 = v[i][j][isSolidIJKplus ? k : k + 1];
				xValue4 = v[i][j + 1][isSolidIJKplus ? k : k + 1];
				yValue1 = v[isSolidIminusJK ? i : i - 1][j][k];
				yValue2 = v[isSolidIminusJK ? i : i - 1][j + 1][k];
				yValue3 = v[isSolidIminusJK ? i : i - 1][j][isSolidIJKplus ? k : k + 1];
				yValue4 = v[isSolidIminusJK ? i : i - 1][j + 1][isSolidIJKplus ? k : k + 1];
			}
			else
			{
				isSolidIJKminus = isSolidCell(i, j, k-1);
				isSolidIminusJK = isSolidCell(i-1, j, k);
				xValue3 = v[i][j][isSolidIJKminus ? k : k - 1];
				xValue4 = v[i][j + 1][isSolidIJKminus ? k : k - 1];
				yValue1 = v[isSolidIminusJK ? i : i - 1][j][k];
				yValue2 = v[isSolidIminusJK ? i : i - 1][j + 1][k];
				yValue3 = v[isSolidIminusJK ? i : i - 1][j][isSolidIJKminus ? k : k - 1];
				yValue4 = v[isSolidIminusJK ? i : i - 1][j + 1][isSolidIJKminus ? k : k - 1];
			}
		} 
		resultV = TrilinearInterpolate<DOUBLE>(xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, yAlpha, zAlpha, xAlpha);
	}
	return resultV;
}

DOUBLE GridField::SampleWWithLocalGrid(const vector3& localPos, INT i, INT j, INT k) const
{
	DOUBLE resultW;
	const vector3 minPt = scnBox.GetMin();
	const vector3 maxPt = scnBox.GetMax();

	INT numX = numXElement;
	INT numY = numYElement;
	INT numZ = numZElement;
	//INT i = 0, j = 0, k = 0;

	//vector3 localPos = fromWorldToLocalGrid(pos, (UINT&)i, (UINT&)j, (UINT&)k);

	DOUBLE xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, xAlpha, yAlpha, zAlpha; 
	
	{
		DOUBLE invertStepSize = 1.0 / stepSize;
		xValue1 = w[i][j][k];
		xValue2 = w[i][j][k + 1];
		zAlpha = localPos.z * invertStepSize;
		xAlpha = signedAbs(0.5 - localPos.x * invertStepSize);
		yAlpha = signedAbs(0.5 - localPos.y * invertStepSize);

		DOUBLE halfStepSize = 0.5 * stepSize;
		DOUBLE isSolidIplusJK, isSolidIminusJK, isSolidIJminusK, isSolidIJplusK;
		if (localPos.y >= halfStepSize)
		{
			if (localPos.x >= halfStepSize)
			{
				isSolidIplusJK = isSolidCell(i+1, j, k);
				isSolidIJplusK = isSolidCell(i, j+1, k);
				xValue3 = w[isSolidIplusJK ? i : i + 1][j][k];
				xValue4 = w[isSolidIplusJK ? i : i + 1][j][k + 1];
				yValue1 = w[i][isSolidIJplusK ? j : j + 1][k];
				yValue2 = w[i][isSolidIJplusK ? j : j + 1][k + 1];
				yValue3 = w[isSolidIplusJK ? i : i + 1][isSolidIJplusK ? j : j + 1][k];
				yValue4 = w[isSolidIplusJK ? i : i + 1][isSolidIJplusK ? j : j + 1][k + 1];
			}
			else
			{
				isSolidIminusJK = isSolidCell(i-1, j, k);
				isSolidIJplusK  = isSolidCell(i, j+1, k);
				xValue3 = w[isSolidIminusJK ? i : i - 1][j][k];
				xValue4 = w[isSolidIminusJK ? i : i - 1][j][k + 1];
				yValue1 = w[i][isSolidIJplusK ? j : j + 1][k];
				yValue2 = w[i][isSolidIJplusK ? j : j + 1][k + 1];
				yValue3 = w[isSolidIminusJK ? i : i - 1][isSolidIJplusK ? j : j + 1][k];
				yValue4 = w[isSolidIminusJK ? i : i - 1][isSolidIJplusK ? j : j + 1][k + 1];
			}
		}
		else 
		{	
			if (localPos.x >= halfStepSize)
			{
				isSolidIplusJK  = isSolidCell(i+1, j, k);
				isSolidIJminusK = isSolidCell(i, j-1, k);
				xValue3 = w[isSolidIplusJK ? i : i + 1][j][k];
				xValue4 = w[isSolidIplusJK ? i : i + 1][j][k + 1];
				yValue1 = w[i][isSolidIJminusK ? j : j - 1][k];
				yValue2 = w[i][isSolidIJminusK ? j : j - 1][k + 1];
				yValue3 = w[isSolidIplusJK ? i : i + 1][isSolidIJminusK ? j : j - 1][k];
				yValue4 = w[isSolidIplusJK ? i : i + 1][isSolidIJminusK ? j : j - 1][k + 1];
			}
			else 
			{
				isSolidIminusJK = isSolidCell(i-1, j, k);
				isSolidIJminusK = isSolidCell(i, j-1, k);
				xValue3 = w[isSolidIminusJK ? i : i - 1][j][k];
				xValue4 = w[isSolidIminusJK ? i : i - 1][j][k + 1];
				yValue1 = w[i][isSolidIJminusK ? j : j - 1][k];
				yValue2 = w[i][isSolidIJminusK ? j : j - 1][k + 1];
				yValue3 = w[isSolidIminusJK ? i : i - 1][isSolidIJminusK ? j : j - 1][k];
				yValue4 = w[isSolidIminusJK ? i : i - 1][isSolidIJminusK ? j : j - 1][k + 1];
			}
		}

		resultW = TrilinearInterpolate<DOUBLE>(xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, zAlpha, xAlpha, yAlpha);
	}
	return resultW;
}


DOUBLE GridField::SampleUWithLocalGrid(const vector3& pos, const vector3& localPos, INT i, INT j, INT k, GridField* coarseLayer, BOOL bPrintf /* = FALSE */) const
{
	if (this == coarseLayer || !coarseLayer)
		return SampleUWithLocalGrid(localPos, i, j, k, bPrintf);
	else
	{
		DOUBLE resultU;
		vector3 localPos;
		INT numX = numXElement;
		INT numY = numYElement;
		INT numZ = numZElement;
		INT i = 0, j = 0, k = 0;

		localPos = fromWorldToLocalGrid(pos, (UINT&)i, (UINT&)j, (UINT&)k);

		DOUBLE xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, xAlpha, yAlpha, zAlpha; 

		
		{
			xValue1 = isVelocityValid('u', i, j, k) ? u[i][j][k] : coarseLayer->SampleU(GetGridCenter(i, j, k) + vector3(-stepSize / 2, 0, 0));
			xValue2 = isVelocityValid('u', i + 1, j, k) ? u[i + 1][j][k] : coarseLayer->SampleU(GetGridCenter(i, j, k) + vector3(stepSize / 2, 0, 0));
			xAlpha = localPos.x / stepSize;
			zAlpha = signedAbs(stepSize / 2 - localPos.z) / stepSize;
			yAlpha = signedAbs(stepSize / 2 - localPos.y) / stepSize;

			if (localPos.y >= 0.5 * stepSize && localPos.z >= 0.5 * stepSize)
			{
				xValue3 = isVelocityValid('u', i, j, k + 1) ? u[i][j][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleU(GetGridCenter(i, j, k + 1) + vector3(-stepSize / 2, 0, 0));
				xValue4 = isVelocityValid('u', i + 1, j, k + 1) ? u[i + 1][j][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleU(GetGridCenter(i, j, k + 1) + vector3(stepSize / 2, 0, 0));
				yValue1 = isVelocityValid('u', i, j + 1, k) ? u[i][isSolidCell(i, j + 1, k) ? j : j + 1][k] : coarseLayer->SampleU(GetGridCenter(i, j + 1, k) + vector3(-stepSize / 2, 0, 0));
				yValue2 = isVelocityValid('u', i + 1, j + 1, k) ? u[i + 1][isSolidCell(i, j + 1, k) ? j : j + 1][k] : coarseLayer->SampleU(GetGridCenter(i, j + 1, k) + vector3(stepSize / 2, 0, 0));
				yValue3 = isVelocityValid('u', i, j + 1, k + 1) ? u[i][isSolidCell(i, j + 1, k) ? j : j + 1][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleU(GetGridCenter(i, j + 1, k + 1) + vector3(-stepSize / 2, 0, 0));
				yValue4 = isVelocityValid('u', i + 1, j + 1, k + 1) ? u[i + 1][isSolidCell(i, j + 1, k) ? j : j + 1][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleU(GetGridCenter(i, j + 1, k + 1) + vector3(stepSize / 2, 0, 0));
			}
			else if (localPos.y < 0.5 * stepSize && localPos.z >= 0.5 * stepSize)
			{
				xValue3 = isVelocityValid('u', i, j, k + 1) ? u[i][j][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleU(GetGridCenter(i, j, k + 1) + vector3(-stepSize / 2, 0, 0));
				xValue4 = isVelocityValid('u', i + 1, j, k + 1) ? u[i + 1][j][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleU(GetGridCenter(i, j, k + 1) + vector3(stepSize / 2, 0, 0));
				yValue1 = isVelocityValid('u', i, j - 1, k) ? u[i][isSolidCell(i, j - 1, k) ? j : j - 1][k] : coarseLayer->SampleU(GetGridCenter(i, j - 1, k) + vector3(-stepSize / 2, 0, 0));
				yValue2 = isVelocityValid('u', i + 1, j - 1, k) ? u[i + 1][isSolidCell(i, j - 1, k) ? j : j - 1][k] : coarseLayer->SampleU(GetGridCenter(i, j - 1, k) + vector3(stepSize / 2, 0, 0));
				yValue3 = isVelocityValid('u', i, j - 1, k + 1) ? u[i][isSolidCell(i, j - 1, k) ? j : j - 1][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleU(GetGridCenter(i, j - 1, k + 1) + vector3(-stepSize / 2, 0, 0));
				yValue4 = isVelocityValid('u', i + 1, j - 1, k + 1) ? u[i + 1][isSolidCell(i, j - 1, k) ? j : j - 1][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleU(GetGridCenter(i, j - 1, k + 1) + vector3(stepSize / 2, 0, 0));
			}
			else if (localPos.y < 0.5 * stepSize && localPos.z < 0.5 * stepSize)
			{
				xValue3 = isVelocityValid('u', i, j, k - 1) ? u[i][j][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleU(GetGridCenter(i, j, k - 1) + vector3(-stepSize / 2, 0, 0));
				xValue4 = isVelocityValid('u', i + 1, j, k - 1) ? u[i + 1][j][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleU(GetGridCenter(i, j, k - 1) + vector3(stepSize / 2, 0, 0));
				yValue1 = isVelocityValid('u', i, j - 1, k) ? u[i][isSolidCell(i, j - 1, k) ? j : j - 1][k] : coarseLayer->SampleU(GetGridCenter(i, j - 1, k) + vector3(-stepSize / 2, 0, 0));
				yValue2 = isVelocityValid('u', i + 1, j - 1, k) ? u[i + 1][isSolidCell(i, j - 1, k) ? j : j - 1][k] : coarseLayer->SampleU(GetGridCenter(i, j - 1, k) + vector3(stepSize / 2, 0, 0));
				yValue3 = isVelocityValid('u', i, j - 1, k - 1) ? u[i][isSolidCell(i, j - 1, k) ? j : j - 1][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleU(GetGridCenter(i, j - 1, k - 1) + vector3(-stepSize / 2, 0, 0));
				yValue4 = isVelocityValid('u', i + 1, j - 1, k - 1) ? u[i + 1][isSolidCell(i, j - 1, k) ? j : j - 1][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleU(GetGridCenter(i, j - 1, k - 1) + vector3(stepSize / 2, 0, 0));
			}
			else if (localPos.y >= 0.5 * stepSize && localPos.z < 0.5 * stepSize)
			{
				xValue3 = isVelocityValid('u', i, j, k - 1) ? u[i][j][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleU(GetGridCenter(i, j, k - 1) + vector3(-stepSize / 2, 0, 0));
				xValue4 = isVelocityValid('u', i + 1, j, k - 1) ? u[i + 1][j][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleU(GetGridCenter(i, j, k - 1) + vector3(stepSize / 2, 0, 0));
				yValue1 = isVelocityValid('u', i, j + 1, k) ? u[i][isSolidCell(i, j + 1, k) ? j : j + 1][k] : coarseLayer->SampleU(GetGridCenter(i, j + 1, k) + vector3(-stepSize / 2, 0, 0));
				yValue2 = isVelocityValid('u', i + 1, j + 1, k) ? u[i + 1][isSolidCell(i, j + 1, k) ? j : j + 1][k] : coarseLayer->SampleU(GetGridCenter(i, j + 1, k) + vector3(stepSize / 2, 0, 0));
				yValue3 = isVelocityValid('u', i, j + 1, k - 1) ? u[i][isSolidCell(i, j + 1, k) ? j : j + 1][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleU(GetGridCenter(i, j + 1, k - 1) + vector3(-stepSize / 2, 0, 0));
				yValue4 = isVelocityValid('u', i + 1, j + 1, k - 1) ? u[i + 1][isSolidCell(i, j + 1, k) ? j : j + 1][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleU(GetGridCenter(i, j + 1, k - 1) + vector3(stepSize / 2, 0, 0));
			}

			resultU = TrilinearInterpolate<DOUBLE>(xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, xAlpha, zAlpha, yAlpha);
		}
		if (bPrintf)
		{
			(*DecoLogger::GetSingleton()) << "xValue1 :" << xValue1 << "\n";
			(*DecoLogger::GetSingleton()) << "xValue2 :" << xValue2 << "\n";
			(*DecoLogger::GetSingleton()) << "xValue3 :" << xValue3 << "\n";
			(*DecoLogger::GetSingleton()) << "xValue4 :" << xValue4 << "\n";
			(*DecoLogger::GetSingleton()) << "yValue1 :" << yValue1 << "\n";
			(*DecoLogger::GetSingleton()) << "yValue2 :" << yValue2 << "\n";
			(*DecoLogger::GetSingleton()) << "yValue3 :" << yValue3 << "\n";
			(*DecoLogger::GetSingleton()) << "yValue4 :" << yValue4 << "\n";
			(*DecoLogger::GetSingleton()) << "xAlpha :" << xAlpha << "\n";
			(*DecoLogger::GetSingleton()) << "yAlpha :" << yAlpha << "\n";
			(*DecoLogger::GetSingleton()) << "zAlpha :" << zAlpha << "\n";
			(*DecoLogger::GetSingleton()) << "resultU :" << resultU << "\n";
		}
		return resultU;
	}
}
DOUBLE GridField::SampleVWithLocalGrid(const vector3& pos, const vector3& localPos, INT i, INT j, INT k, GridField* coarseLayer) const
{
	if (this == coarseLayer || !coarseLayer)
		return SampleVWithLocalGrid(localPos, i, j, k);
	else
	{
		DOUBLE resultV;
		vector3 localPos;
		INT numX = numXElement;
		INT numY = numYElement;
		INT numZ = numZElement;
		INT i = 0, j = 0, k = 0;

		localPos = fromWorldToLocalGrid(pos, (UINT&)i, (UINT&)j, (UINT&)k);

		DOUBLE xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, xAlpha, yAlpha, zAlpha; 

		{
			xValue1 = isVelocityValid('v', i, j, k) ? v[i][j][k] : coarseLayer->SampleV(GetGridCenter(i, j, k) + vector3(0, -stepSize / 2, 0));
			xValue2 = isVelocityValid('v', i, j + 1, k) ? v[i][j + 1][k] : coarseLayer->SampleV(GetGridCenter(i, j, k) + vector3(0, stepSize / 2, 0));
			yAlpha = localPos.y / stepSize;
			zAlpha = signedAbs(stepSize / 2 - localPos.z) / stepSize;
			xAlpha = signedAbs(stepSize / 2 - localPos.x) / stepSize;

			if (localPos.x >= 0.5 * stepSize && localPos.z >= 0.5 * stepSize)
			{
				xValue3 = (isVelocityValid('v', i, j, k + 1) ? v[i][j][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleV(GetGridCenter(i, j, k + 1) + vector3(0, -stepSize / 2, 0)));
				xValue4 = (isVelocityValid('v', i, j + 1, k + 1) ? v[i][j + 1][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleV(GetGridCenter(i, j, k + 1) + vector3(0, stepSize / 2, 0)));
				yValue1 = (isVelocityValid('v', i + 1, j, k) ? v[isSolidCell(i + 1, j, k) ? i : i + 1][j][k] : coarseLayer->SampleV(GetGridCenter(i + 1, j, k) + vector3(0, -stepSize / 2, 0)));
				yValue2 = (isVelocityValid('v', i + 1, j + 1, k) ? v[isSolidCell(i + 1, j, k) ? i : i + 1][j + 1][k] : coarseLayer->SampleV(GetGridCenter(i + 1, j, k) + vector3(0, stepSize / 2, 0)));
				yValue3 = (isVelocityValid('v', i + 1, j, k + 1) ? v[isSolidCell(i + 1, j, k) ? i : i + 1][j][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleV(GetGridCenter(i + 1, j, k + 1) + vector3(0, -stepSize / 2,0)));
				yValue4 = (isVelocityValid('v', i + 1, j + 1, k + 1) ? v[isSolidCell(i + 1, j, k) ? i : i + 1][j + 1][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleV(GetGridCenter(i + 1, j , k + 1) + vector3(0, stepSize / 2, 0)));
			}
			else if (localPos.x < 0.5 * stepSize && localPos.z >= 0.5 * stepSize)
			{
				xValue3 = (isVelocityValid('v', i, j, k + 1) ? v[i][j][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleV(GetGridCenter(i, j, k + 1) + vector3(0, -stepSize / 2, 0)));
				xValue4 = (isVelocityValid('v', i, j + 1, k + 1) ? v[i][j + 1][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleV(GetGridCenter(i, j, k + 1) + vector3(0, stepSize / 2, 0)));
				yValue1 = (isVelocityValid('v', i - 1, j, k) ? v[isSolidCell(i - 1, j, k) ? i : i - 1][j][k] : coarseLayer->SampleV(GetGridCenter(i - 1, j, k) + vector3(0, -stepSize / 2, 0)));
				yValue2 = (isVelocityValid('v', i - 1, j + 1, k) ? v[isSolidCell(i - 1, j, k) ? i : i - 1][j + 1][k] : coarseLayer->SampleV(GetGridCenter(i - 1, j, k) + vector3(0, stepSize / 2, 0)));
				yValue3 = (isVelocityValid('v', i - 1, j, k + 1) ? v[isSolidCell(i - 1, j, k) ? i : i - 1][j][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleV(GetGridCenter(i - 1, j, k + 1) + vector3(0, -stepSize / 2, 0)));
				yValue4 = (isVelocityValid('v', i - 1, j + 1, k + 1) ? v[isSolidCell(i - 1, j, k) ? i : i - 1][j + 1][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleV(GetGridCenter(i - 1, j, k + 1) + vector3(0, stepSize / 2, 0)));
			}
			else if (localPos.x < 0.5 * stepSize && localPos.z < 0.5 * stepSize)
			{
				xValue3 = (isVelocityValid('v', i, j, k - 1) ? v[i][j][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleV(GetGridCenter(i, j, k - 1) + vector3(0, -stepSize / 2, 0)));
				xValue4 = (isVelocityValid('v', i, j + 1, k - 1) ? v[i][j + 1][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleV(GetGridCenter(i, j, k - 1) + vector3(0, stepSize / 2, 0)));
				yValue1 = (isVelocityValid('v', i - 1, j, k) ? v[isSolidCell(i - 1, j, k) ? i : i - 1][j][k] : coarseLayer->SampleV(GetGridCenter(i - 1, j, k) + vector3(0, -stepSize / 2, 0)));
				yValue2 = (isVelocityValid('v', i - 1, j + 1, k) ? v[isSolidCell(i - 1, j, k) ? i : i - 1][j + 1][k] : coarseLayer->SampleV(GetGridCenter(i - 1, j, k) + vector3(0, stepSize / 2, 0)));
				yValue3 = (isVelocityValid('v', i - 1, j, k - 1) ? v[isSolidCell(i - 1, j, k) ? i : i - 1][j][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleV(GetGridCenter(i - 1, j, k - 1) + vector3(0, -stepSize / 2, 0))); 
				yValue4 = (isVelocityValid('v', i - 1, j + 1, k - 1) ? v[isSolidCell(i - 1, j, k) ? i : i - 1][j + 1][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleV(GetGridCenter(i - 1, j, k - 1) + vector3(0, stepSize / 2, 0)));
			}
			else if (localPos.x >= 0.5 * stepSize && localPos.z < 0.5 * stepSize)
			{
				xValue3 = (isVelocityValid('v', i, j, k - 1) ? v[i][j][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleV(GetGridCenter(i, j, k - 1) + vector3(0, -stepSize / 2, 0)));
				xValue4 = (isVelocityValid('v', i, j + 1, k - 1) ? v[i][j + 1][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleV(GetGridCenter(i, j, k - 1) + vector3(0, stepSize / 2, 0)));
				yValue1 = (isVelocityValid('v', i + 1, j, k) ? v[isSolidCell(i + 1, j, k) ? i : i + 1][j][k] : coarseLayer->SampleV(GetGridCenter(i + 1, j, k) + vector3(0, -stepSize / 2, 0)));
				yValue2 = (isVelocityValid('v', i + 1, j + 1, k) ? v[isSolidCell(i + 1, j, k) ? i : i + 1][j + 1][k] : coarseLayer->SampleV(GetGridCenter(i + 1, j, k) + vector3(0, stepSize / 2, 0)));
				yValue3 = (isVelocityValid('v', i + 1, j, k - 1) ? v[isSolidCell(i + 1, j, k) ? i : i + 1][j][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleV(GetGridCenter(i + 1, j, k - 1) + vector3(0, - stepSize / 2, 0)));
				yValue4 = (isVelocityValid('v', i + 1, j + 1, k - 1) ? v[isSolidCell(i + 1, j, k) ? i : i + 1][j + 1][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleV(GetGridCenter(i + 1, j, k - 1) + vector3(0, stepSize / 2, 0)));
			}
			else 
				assert(0);
			resultV = TrilinearInterpolate<DOUBLE>(xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, yAlpha, zAlpha, xAlpha);
		}
		return resultV;

	}
}
DOUBLE GridField::SampleWWithLocalGrid(const vector3& pos, const vector3& localPos, INT i, INT j, INT k, GridField* coarseLayer) const
{
	if (this == coarseLayer || !coarseLayer)
		return SampleWWithLocalGrid(localPos, i, j, k);
	else
	{
		DOUBLE resultW;
		vector3 localPos;
		INT numX = numXElement;
		INT numY = numYElement;
		INT numZ = numZElement;
		INT i = 0, j = 0, k = 0;

		localPos = fromWorldToLocalGrid(pos, (UINT&)i, (UINT&)j, (UINT&)k);

		DOUBLE xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, xAlpha, yAlpha, zAlpha; 
		
		{
			xValue1 = isVelocityValid('w', i, j, k) ? w[i][j][k] : coarseLayer->SampleW(GetGridCenter(i, j, k) + vector3(0, 0, -stepSize / 2));
			xValue2 = isVelocityValid('w', i, j, k + 1) ? w[i][j][k + 1] : coarseLayer->SampleW(GetGridCenter(i, j, k) + vector3(0, 0, stepSize / 2));
			zAlpha = localPos.z / stepSize;
			xAlpha = signedAbs(stepSize / 2 - localPos.x) / stepSize;
			yAlpha = signedAbs(stepSize / 2 - localPos.y) / stepSize;

			if (localPos.y >= 0.5 * stepSize && localPos.x >= 0.5 * stepSize)
			{
				xValue3 = isVelocityValid('w', i + 1, j, k) ? w[isSolidCell(i + 1, j, k) ? i : i + 1][j][k] : coarseLayer->SampleW(GetGridCenter(i + 1,j ,k) + vector3(0, 0, -stepSize / 2));
				xValue4 = isVelocityValid('w', i + 1, j, k + 1) ? w[isSolidCell(i + 1, j, k) ? i : i + 1][j][k + 1] : coarseLayer->SampleW(GetGridCenter(i + 1, j, k) + vector3(0, 0, stepSize / 2));
				yValue1 = isVelocityValid('w', i, j + 1, k) ? w[i][isSolidCell(i, j + 1, k) ? j : j + 1][k] : coarseLayer->SampleW(GetGridCenter(i, j + 1, k) + vector3(0, 0, -stepSize / 2));
				yValue2 = isVelocityValid('w', i, j + 1, k + 1) ? w[i][isSolidCell(i, j + 1, k) ? j : j + 1][k + 1] : coarseLayer->SampleW(GetGridCenter(i, j + 1, k) + vector3(0, 0, stepSize / 2));
				yValue3 = isVelocityValid('w', i + 1, j + 1, k) ? w[isSolidCell(i + 1, j, k) ? i : i + 1][isSolidCell(i, j + 1, k) ? j : j + 1][k] : coarseLayer->SampleW(GetGridCenter(i + 1, j + 1, k) + vector3(0, 0, -stepSize / 2));
				yValue4 = isVelocityValid('w', i + 1, j + 1, k + 1) ? w[isSolidCell(i + 1, j, k) ? i : i + 1][isSolidCell(i, j + 1, k) ? j : j + 1][k + 1] : coarseLayer->SampleW(GetGridCenter(i + 1, j + 1, k) + vector3(0, 0, stepSize / 2));
			}
			else if (localPos.y < 0.5 * stepSize && localPos.x >= 0.5 * stepSize)
			{
				xValue3 = isVelocityValid('w', i + 1, j, k) ? w[isSolidCell(i + 1, j, k) ? i : i + 1][j][k] : coarseLayer->SampleW(GetGridCenter(i + 1,j ,k) + vector3(0, 0, -stepSize / 2));
				xValue4 = isVelocityValid('w', i + 1, j, k + 1) ? w[isSolidCell(i + 1, j, k) ? i : i + 1][j][k + 1] : coarseLayer->SampleW(GetGridCenter(i + 1, j, k) + vector3(0, 0, stepSize / 2));
				yValue1 = isVelocityValid('w', i, j - 1, k) ? w[i][isSolidCell(i, j - 1, k) ? j : j - 1][k] : coarseLayer->SampleW(GetGridCenter(i, j - 1, k) + vector3(0, 0, -stepSize / 2));
				yValue2 = isVelocityValid('w', i, j - 1, k + 1) ? w[i][isSolidCell(i, j - 1, k) ? j : j - 1][k + 1] : coarseLayer->SampleW(GetGridCenter(i, j - 1, k) + vector3(0, 0, stepSize / 2));
				yValue3 = isVelocityValid('w', i + 1, j - 1, k) ? w[isSolidCell(i + 1, j, k) ? i : i + 1][isSolidCell(i, j - 1, k) ? j : j - 1][k] : coarseLayer->SampleW(GetGridCenter(i + 1, j - 1, k) + vector3(0, 0, -stepSize / 2));
				yValue4 = isVelocityValid('w', i + 1, j - 1, k + 1) ? w[isSolidCell(i + 1, j, k) ? i : i + 1][isSolidCell(i, j - 1, k) ? j : j - 1][k + 1] : coarseLayer->SampleW(GetGridCenter(i + 1, j - 1, k) + vector3(0, 0, stepSize / 2));
			}
			else if (localPos.y < 0.5 * stepSize && localPos.x < 0.5 * stepSize)
			{
				xValue3 = isVelocityValid('w', i - 1, j, k) ? w[isSolidCell(i - 1, j, k) ? i : i - 1][j][k] : coarseLayer->SampleW(GetGridCenter(i - 1, j, k) + vector3(0, 0, -stepSize / 2));
				xValue4 = isVelocityValid('w', i - 1, j, k + 1) ? w[isSolidCell(i - 1, j, k) ? i : i - 1][j][k + 1] : coarseLayer->SampleW(GetGridCenter(i - 1, j, k) + vector3(0, 0, stepSize / 2));
				yValue1 = isVelocityValid('w', i, j - 1, k) ? w[i][isSolidCell(i, j - 1, k) ? j : j - 1][k] : coarseLayer->SampleW(GetGridCenter(i, j - 1, k) + vector3(0, 0, -stepSize / 2));
				yValue2 = isVelocityValid('w', i, j - 1, k + 1) ? w[i][isSolidCell(i, j - 1, k) ? j : j - 1][k + 1] : coarseLayer->SampleW(GetGridCenter(i, j - 1, k) + vector3(0, 0, stepSize / 2));
				yValue3 = isVelocityValid('w', i - 1, j - 1, k) ? w[isSolidCell(i - 1, j, k) ? i : i - 1][isSolidCell(i, j - 1, k) ? j : j - 1][k] : coarseLayer->SampleW(GetGridCenter(i - 1, j - 1, k) + vector3(0, 0, -stepSize / 2));
				yValue4 = isVelocityValid('w', i - 1, j - 1, k + 1) ? w[isSolidCell(i - 1, j, k) ? i : i - 1][isSolidCell(i, j - 1, k) ? j : j - 1][k + 1] : coarseLayer->SampleW(GetGridCenter(i - 1, j - 1, k) + vector3(0, 0, stepSize / 2));
			}
			else if (localPos.y >= 0.5 * stepSize && localPos.x < 0.5 * stepSize)
			{
				xValue3 = isVelocityValid('w', i - 1, j, k) ? w[isSolidCell(i - 1, j, k) ? i : i - 1][j][k] : coarseLayer->SampleW(GetGridCenter(i - 1, j, k) + vector3(0, 0, -stepSize / 2));
				xValue4 = isVelocityValid('w', i - 1, j, k + 1) ? w[isSolidCell(i - 1, j, k) ? i : i - 1][j][k + 1] : coarseLayer->SampleW(GetGridCenter(i - 1, j, k) + vector3(0, 0, stepSize / 2));
				yValue1 = isVelocityValid('w', i, j + 1, k) ? w[i][isSolidCell(i, j + 1, k) ? j : j + 1][k] : coarseLayer->SampleW(GetGridCenter(i, j + 1, k) + vector3(0, 0, -stepSize / 2));
				yValue2 = isVelocityValid('w', i, j + 1, k + 1) ? w[i][isSolidCell(i, j + 1, k) ? j : j + 1][k + 1] : coarseLayer->SampleW(GetGridCenter(i, j + 1, k) + vector3(0, 0, stepSize / 2));
				yValue3 = isVelocityValid('w', i - 1, j + 1, k) ? w[isSolidCell(i - 1, j, k) ? i : i - 1][isSolidCell(i, j + 1, k) ? j : j + 1][k] : coarseLayer->SampleW(GetGridCenter(i - 1, j + 1, k) + vector3(0, 0, -stepSize / 2));
				yValue4 = isVelocityValid('w', i - 1, j + 1, k + 1) ? w[isSolidCell(i - 1, j, k) ? i : i - 1][isSolidCell(i, j + 1, k) ? j : j + 1][k + 1] : coarseLayer->SampleW(GetGridCenter(i - 1, j + 1, k) + vector3(0, 0, stepSize / 2));
			}
			else
				assert(0);
			resultW = TrilinearInterpolate<DOUBLE>(xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, zAlpha, xAlpha, yAlpha);
		}
		return resultW;
	}
}

//  [4/24/2008 HATEVOL] ADDED end

DOUBLE GridField::SampleU(const vector3& pos, BOOL bPrint) const
{
	DOUBLE resultU;
	const vector3 minPt = scnBox.GetMin();
	const vector3 maxPt = scnBox.GetMax();

	
	INT numX = numXElement;
	INT numY = numYElement;
	INT numZ = numZElement;
	INT i = 0, j = 0, k = 0;

	vector3 localPos = fromWorldToLocalGrid(pos, (UINT&)i, (UINT&)j, (UINT&)k);

	DOUBLE xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, xAlpha, yAlpha, zAlpha; 


	{
		DOUBLE invertStepSize = 1 / stepSize;
		xValue1 = u[i][j][k];
		xValue2 = u[i + 1][j][k];
		xAlpha = localPos.x * invertStepSize;
		zAlpha = signedAbs(0.5 - localPos.z * invertStepSize );
		yAlpha = signedAbs(0.5 - localPos.y * invertStepSize);


		DOUBLE halfStepSize = 0.5 * stepSize;
		BOOL isSolidIJKminus, isSolidIJKplus, isSolidIJminusK, isSolidIJplusK;
		
		if ( localPos.y >= halfStepSize )
		{
		    if ( localPos.z >= halfStepSize )
			{
			    isSolidIJKplus  = isSolidCell(i, j, k+1);
				isSolidIJplusK  = isSolidCell(i, j+1, k);

				xValue3 = u[i][j][isSolidIJKplus ? k : k + 1];
				yValue1 = u[i][isSolidIJplusK ? j : j + 1][k];
				yValue3 = u[i][isSolidIJplusK ? j : j + 1][isSolidIJKplus ? k : k + 1];
				xValue4 = u[i + 1][j][isSolidIJKplus ? k : k + 1];
				yValue2 = u[i + 1][isSolidIJplusK ? j : j + 1][k];
				yValue4 = u[i + 1][isSolidIJplusK ? j : j + 1][isSolidIJKplus ? k : k + 1];

			}
			else {
				isSolidIJKminus = isSolidCell(i, j, k-1);
				isSolidIJplusK  = isSolidCell(i, j+1, k);

				xValue3 = u[i][j][isSolidIJKminus ? k : k - 1];
				yValue1 = u[i][isSolidIJplusK ? j : j + 1][k];
				yValue3 = u[i][isSolidIJplusK ? j : j + 1][isSolidIJKminus ? k : k - 1];
				xValue4 = u[i + 1][j][isSolidIJKminus ? k : k - 1];
				yValue2 = u[i + 1][isSolidIJplusK ? j : j + 1][k];
				yValue4 = u[i + 1][isSolidIJplusK ? j : j + 1][isSolidIJKminus ? k : k - 1];
			}
		}
		else 
		{
			if ( localPos.z >= halfStepSize)
		    {
				isSolidIJKplus  = isSolidCell(i, j, k+1);
				isSolidIJminusK = isSolidCell(i, j-1, k);

				xValue3 = u[i][j][isSolidIJKplus ? k : k + 1];
				yValue1 = u[i][isSolidIJminusK ? j : j - 1][k];
				yValue3 = u[i][isSolidIJminusK ? j : j - 1][isSolidIJKplus ? k : k + 1];
				xValue4 = u[i + 1][j][isSolidIJKplus ? k : k + 1]; 
				yValue2 = u[i + 1][isSolidIJminusK ? j : j - 1][k];
				yValue4 = u[i + 1][isSolidIJminusK ? j : j - 1][isSolidIJKplus ? k : k + 1];
			}
		    else 
		    {
				isSolidIJKminus = isSolidCell(i, j, k-1);
				isSolidIJminusK = isSolidCell(i, j-1, k);
		    
				xValue3 = u[i][j][isSolidIJKminus ? k : k - 1];
				yValue1 = u[i][isSolidIJminusK ? j : j - 1][k];
				yValue3 = u[i][isSolidIJminusK ? j : j - 1][isSolidIJKminus ? k : k - 1];
				xValue4 = u[i + 1][j][isSolidIJKminus ? k : k - 1];
				yValue2 = u[i + 1][isSolidIJminusK ? j : j - 1][k];
				yValue4 = u[i + 1][isSolidIJminusK ? j : j - 1][isSolidIJKminus ? k : k - 1];

			}
			
		}
		//  [4/8/2008 HATEVOL] ADDED end
		//else 
		resultU = TrilinearInterpolate<DOUBLE>(xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, xAlpha, zAlpha, yAlpha);
	}
	if (bPrint)
	{
		(*DecoLogger::GetSingleton()) << "xValue1 :" << xValue1 << "\n";
		(*DecoLogger::GetSingleton()) << "xValue2 :" << xValue2 << "\n";
		(*DecoLogger::GetSingleton()) << "xValue3 :" << xValue3 << "\n";
		(*DecoLogger::GetSingleton()) << "xValue4 :" << xValue4 << "\n";
		(*DecoLogger::GetSingleton()) << "yValue1 :" << yValue1 << "\n";
		(*DecoLogger::GetSingleton()) << "yValue2 :" << yValue2 << "\n";
		(*DecoLogger::GetSingleton()) << "yValue3 :" << yValue3 << "\n";
		(*DecoLogger::GetSingleton()) << "yValue4 :" << yValue4 << "\n";
		(*DecoLogger::GetSingleton()) << "xAlpha :" << xAlpha << "\n";
		(*DecoLogger::GetSingleton()) << "yAlpha :" << yAlpha << "\n";
		(*DecoLogger::GetSingleton()) << "zAlpha :" << zAlpha << "\n";
		(*DecoLogger::GetSingleton()) << "resultU :" << resultU << "\n";
	}
	return resultU;
}

DOUBLE GridField::SampleV(const vector3& pos) const
{
	DOUBLE resultV;
	const vector3 minPt = scnBox.GetMin();
	const vector3 maxPt = scnBox.GetMax();

	
	INT numX = numXElement;
	INT numY = numYElement;
	INT numZ = numZElement;
	INT i = 0, j = 0, k = 0;

	vector3 localPos = fromWorldToLocalGrid(pos, (UINT&)i, (UINT&)j, (UINT&)k);

	DOUBLE xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, xAlpha, yAlpha, zAlpha; 


	{
		DOUBLE invertStepSize = 1.0 / stepSize;
		xValue1 = v[i][j][k];
		xValue2 = v[i][j + 1][k];
		yAlpha = localPos.y * invertStepSize;
		zAlpha = signedAbs(0.5 - localPos.z * invertStepSize);
		xAlpha = signedAbs(0.5 - localPos.x * invertStepSize);

		DOUBLE halfStepSize = 0.5 * stepSize;
		BOOL isSolidIJKminus, isSolidIJKplus, isSolidIminusJK, isSolidIplusJK;
		if (localPos.x >= halfStepSize)
		{
			if ( localPos.z >= halfStepSize)
		    {
				isSolidIJKplus = isSolidCell(i, j, k+1);
				isSolidIplusJK = isSolidCell(i+1, j, k);
			    xValue3 = v[i][j][isSolidIJKplus ? k : k + 1];
			    xValue4 = v[i][j + 1][isSolidIJKplus ? k : k + 1];
			    yValue1 = v[isSolidIplusJK ? i : i + 1][j][k];
			    yValue2 = v[isSolidIplusJK ? i : i + 1][j + 1][k];
			    yValue3 = v[isSolidIplusJK ? i : i + 1][j][isSolidIJKplus ? k : k + 1];
			    yValue4 = v[isSolidIplusJK ? i : i + 1][j + 1][isSolidIJKplus ? k : k + 1];
		    }
			else
			{
				isSolidIJKminus = isSolidCell(i, j, k-1);
				isSolidIplusJK  = isSolidCell(i+1, j, k);
				xValue3 = v[i][j][isSolidIJKminus ? k : k - 1];
				xValue4 = v[i][j + 1][isSolidIJKminus ? k : k - 1];
				yValue1 = v[isSolidIplusJK ? i : i + 1][j][k];
				yValue2 = v[isSolidIplusJK ? i : i + 1][j + 1][k];
				yValue3 = v[isSolidIplusJK ? i : i + 1][j][isSolidIJKminus ? k : k - 1];
				yValue4 = v[isSolidIplusJK ? i : i + 1][j + 1][isSolidIJKminus ? k : k - 1];
			}
		}
		else 
		{
			if (localPos.z >= halfStepSize)
		    {
				isSolidIJKplus  = isSolidCell(i, j, k+1);
				isSolidIminusJK = isSolidCell(i-1, j, k);
			    xValue3 = v[i][j][isSolidIJKplus ? k : k + 1];
			    xValue4 = v[i][j + 1][isSolidIJKplus ? k : k + 1];
			    yValue1 = v[isSolidIminusJK ? i : i - 1][j][k];
			    yValue2 = v[isSolidIminusJK ? i : i - 1][j + 1][k];
			    yValue3 = v[isSolidIminusJK ? i : i - 1][j][isSolidIJKplus ? k : k + 1];
			    yValue4 = v[isSolidIminusJK ? i : i - 1][j + 1][isSolidIJKplus ? k : k + 1];
		    }
		    else
			{
				isSolidIJKminus = isSolidCell(i, j, k-1);
				isSolidIminusJK = isSolidCell(i-1, j, k);
			    xValue3 = v[i][j][isSolidIJKminus ? k : k - 1];
			    xValue4 = v[i][j + 1][isSolidIJKminus ? k : k - 1];
			    yValue1 = v[isSolidIminusJK ? i : i - 1][j][k];
			    yValue2 = v[isSolidIminusJK ? i : i - 1][j + 1][k];
			    yValue3 = v[isSolidIminusJK ? i : i - 1][j][isSolidIJKminus ? k : k - 1];
			    yValue4 = v[isSolidIminusJK ? i : i - 1][j + 1][isSolidIJKminus ? k : k - 1];
		    }
		} 
		resultV = TrilinearInterpolate<DOUBLE>(xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, yAlpha, zAlpha, xAlpha);
	}
	return resultV;
}

DOUBLE GridField::SampleW(const vector3& pos) const
{
	DOUBLE resultW;
	const vector3 minPt = scnBox.GetMin();
	const vector3 maxPt = scnBox.GetMax();

	
	INT numX = numXElement;
	INT numY = numYElement;
	INT numZ = numZElement;
	INT i = 0, j = 0, k = 0;

	vector3 localPos = fromWorldToLocalGrid(pos, (UINT&)i, (UINT&)j, (UINT&)k);

	DOUBLE xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, xAlpha, yAlpha, zAlpha; 

	{
		DOUBLE invertStepSize = 1.0 / stepSize;
		xValue1 = w[i][j][k];
		xValue2 = w[i][j][k + 1];
		zAlpha = localPos.z * invertStepSize;
		xAlpha = signedAbs(0.5 - localPos.x * invertStepSize);
		yAlpha = signedAbs(0.5 - localPos.y * invertStepSize);
		DOUBLE halfStepSize = 0.5 * stepSize;
		DOUBLE isSolidIplusJK, isSolidIminusJK, isSolidIJminusK, isSolidIJplusK;
		if (localPos.y >= halfStepSize)
		{
			if (localPos.x >= halfStepSize)
		    {
			    isSolidIplusJK = isSolidCell(i+1, j, k);
			    isSolidIJplusK = isSolidCell(i, j+1, k);
			    xValue3 = w[isSolidIplusJK ? i : i + 1][j][k];
			    xValue4 = w[isSolidIplusJK ? i : i + 1][j][k + 1];
			    yValue1 = w[i][isSolidIJplusK ? j : j + 1][k];
			    yValue2 = w[i][isSolidIJplusK ? j : j + 1][k + 1];
			    yValue3 = w[isSolidIplusJK ? i : i + 1][isSolidIJplusK ? j : j + 1][k];
			    yValue4 = w[isSolidIplusJK ? i : i + 1][isSolidIJplusK ? j : j + 1][k + 1];
		    }
			else
			{
				isSolidIminusJK = isSolidCell(i-1, j, k);
				isSolidIJplusK  = isSolidCell(i, j+1, k);
				xValue3 = w[isSolidIminusJK ? i : i - 1][j][k];
				xValue4 = w[isSolidIminusJK ? i : i - 1][j][k + 1];
				yValue1 = w[i][isSolidIJplusK ? j : j + 1][k];
				yValue2 = w[i][isSolidIJplusK ? j : j + 1][k + 1];
				yValue3 = w[isSolidIminusJK ? i : i - 1][isSolidIJplusK ? j : j + 1][k];
				yValue4 = w[isSolidIminusJK ? i : i - 1][isSolidIJplusK ? j : j + 1][k + 1];
			}
		}
		else 
		{	
		    if (localPos.x >= halfStepSize)
		    {
			    isSolidIplusJK  = isSolidCell(i+1, j, k);
			    isSolidIJminusK = isSolidCell(i, j-1, k);
			    xValue3 = w[isSolidIplusJK ? i : i + 1][j][k];
			    xValue4 = w[isSolidIplusJK ? i : i + 1][j][k + 1];
			    yValue1 = w[i][isSolidIJminusK ? j : j - 1][k];
			    yValue2 = w[i][isSolidIJminusK ? j : j - 1][k + 1];
			    yValue3 = w[isSolidIplusJK ? i : i + 1][isSolidIJminusK ? j : j - 1][k];
			    yValue4 = w[isSolidIplusJK ? i : i + 1][isSolidIJminusK ? j : j - 1][k + 1];
		    }
		    else 
		    {
			    isSolidIminusJK = isSolidCell(i-1, j, k);
			    isSolidIJminusK = isSolidCell(i, j-1, k);
			    xValue3 = w[isSolidIminusJK ? i : i - 1][j][k];
			    xValue4 = w[isSolidIminusJK ? i : i - 1][j][k + 1];
			    yValue1 = w[i][isSolidIJminusK ? j : j - 1][k];
			    yValue2 = w[i][isSolidIJminusK ? j : j - 1][k + 1];
			    yValue3 = w[isSolidIminusJK ? i : i - 1][isSolidIJminusK ? j : j - 1][k];
			    yValue4 = w[isSolidIminusJK ? i : i - 1][isSolidIJminusK ? j : j - 1][k + 1];
		    }
		}
		resultW = TrilinearInterpolate<DOUBLE>(xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, zAlpha, xAlpha, yAlpha);
	}
	return resultW;
}

BOOL GridField::isVelocityValid(char dim, INT i, INT j, INT k) const
{
	switch (dim)
	{

	case 'u':
		return (/*isSolidCell(i - 1, j, k) && isSolidCell(i, j, k) || */isInDomain(i - 1, j, k) || isInDomain(i, j, k));
		break;
	case 'v':
		return (/*isSolidCell(i, j - 1, k) && isSolidCell(i, j, k) || */isInDomain(i, j - 1, k) || isInDomain(i, j, k));
		break;
	case 'w':
		return (/*isSolidCell(i, j, k - 1) && isSolidCell(i, j, k) || */isInDomain(i, j, k - 1) || isInDomain(i, j, k));
		break;

	default:
		return FALSE;
	}
}

DOUBLE GridField::SampleU(const vector3& pos, GridField* coarseLayer, BOOL bPrint) const
{
	if (this == coarseLayer || !coarseLayer)
		return SampleU(pos, bPrint);
	else
	{
		DOUBLE resultU;
		vector3 localPos;
		INT numX = numXElement;
		INT numY = numYElement;
		INT numZ = numZElement;
		INT i = 0, j = 0, k = 0;

		localPos = fromWorldToLocalGrid(pos, (UINT&)i, (UINT&)j, (UINT&)k);

		DOUBLE xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, xAlpha, yAlpha, zAlpha; 


		{
			xValue1 = isVelocityValid('u', i, j, k) ? u[i][j][k] : coarseLayer->SampleU(GetGridCenter(i, j, k) + vector3(-stepSize / 2, 0, 0));
			xValue2 = isVelocityValid('u', i + 1, j, k) ? u[i + 1][j][k] : coarseLayer->SampleU(GetGridCenter(i, j, k) + vector3(stepSize / 2, 0, 0));
			xAlpha = localPos.x / stepSize;
			zAlpha = signedAbs(stepSize / 2 - localPos.z) / stepSize;
			yAlpha = signedAbs(stepSize / 2 - localPos.y) / stepSize;

			if (localPos.y >= 0.5 * stepSize && localPos.z >= 0.5 * stepSize)
			{
				xValue3 = isVelocityValid('u', i, j, k + 1) ? u[i][j][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleU(GetGridCenter(i, j, k + 1) + vector3(-stepSize / 2, 0, 0));
				xValue4 = isVelocityValid('u', i + 1, j, k + 1) ? u[i + 1][j][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleU(GetGridCenter(i, j, k + 1) + vector3(stepSize / 2, 0, 0));
				yValue1 = isVelocityValid('u', i, j + 1, k) ? u[i][isSolidCell(i, j + 1, k) ? j : j + 1][k] : coarseLayer->SampleU(GetGridCenter(i, j + 1, k) + vector3(-stepSize / 2, 0, 0));
				yValue2 = isVelocityValid('u', i + 1, j + 1, k) ? u[i + 1][isSolidCell(i, j + 1, k) ? j : j + 1][k] : coarseLayer->SampleU(GetGridCenter(i, j + 1, k) + vector3(stepSize / 2, 0, 0));
				yValue3 = isVelocityValid('u', i, j + 1, k + 1) ? u[i][isSolidCell(i, j + 1, k) ? j : j + 1][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleU(GetGridCenter(i, j + 1, k + 1) + vector3(-stepSize / 2, 0, 0));
				yValue4 = isVelocityValid('u', i + 1, j + 1, k + 1) ? u[i + 1][isSolidCell(i, j + 1, k) ? j : j + 1][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleU(GetGridCenter(i, j + 1, k + 1) + vector3(stepSize / 2, 0, 0));
			}
			else if (localPos.y < 0.5 * stepSize && localPos.z >= 0.5 * stepSize)
			{
				xValue3 = isVelocityValid('u', i, j, k + 1) ? u[i][j][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleU(GetGridCenter(i, j, k + 1) + vector3(-stepSize / 2, 0, 0));
				xValue4 = isVelocityValid('u', i + 1, j, k + 1) ? u[i + 1][j][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleU(GetGridCenter(i, j, k + 1) + vector3(stepSize / 2, 0, 0));
				yValue1 = isVelocityValid('u', i, j - 1, k) ? u[i][isSolidCell(i, j - 1, k) ? j : j - 1][k] : coarseLayer->SampleU(GetGridCenter(i, j - 1, k) + vector3(-stepSize / 2, 0, 0));
				yValue2 = isVelocityValid('u', i + 1, j - 1, k) ? u[i + 1][isSolidCell(i, j - 1, k) ? j : j - 1][k] : coarseLayer->SampleU(GetGridCenter(i, j - 1, k) + vector3(stepSize / 2, 0, 0));
				yValue3 = isVelocityValid('u', i, j - 1, k + 1) ? u[i][isSolidCell(i, j - 1, k) ? j : j - 1][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleU(GetGridCenter(i, j - 1, k + 1) + vector3(-stepSize / 2, 0, 0));
				yValue4 = isVelocityValid('u', i + 1, j - 1, k + 1) ? u[i + 1][isSolidCell(i, j - 1, k) ? j : j - 1][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleU(GetGridCenter(i, j - 1, k + 1) + vector3(stepSize / 2, 0, 0));
			}
			else if (localPos.y < 0.5 * stepSize && localPos.z < 0.5 * stepSize)
			{
				xValue3 = isVelocityValid('u', i, j, k - 1) ? u[i][j][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleU(GetGridCenter(i, j, k - 1) + vector3(-stepSize / 2, 0, 0));
				xValue4 = isVelocityValid('u', i + 1, j, k - 1) ? u[i + 1][j][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleU(GetGridCenter(i, j, k - 1) + vector3(stepSize / 2, 0, 0));
				yValue1 = isVelocityValid('u', i, j - 1, k) ? u[i][isSolidCell(i, j - 1, k) ? j : j - 1][k] : coarseLayer->SampleU(GetGridCenter(i, j - 1, k) + vector3(-stepSize / 2, 0, 0));
				yValue2 = isVelocityValid('u', i + 1, j - 1, k) ? u[i + 1][isSolidCell(i, j - 1, k) ? j : j - 1][k] : coarseLayer->SampleU(GetGridCenter(i, j - 1, k) + vector3(stepSize / 2, 0, 0));
				yValue3 = isVelocityValid('u', i, j - 1, k - 1) ? u[i][isSolidCell(i, j - 1, k) ? j : j - 1][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleU(GetGridCenter(i, j - 1, k - 1) + vector3(-stepSize / 2, 0, 0));
				yValue4 = isVelocityValid('u', i + 1, j - 1, k - 1) ? u[i + 1][isSolidCell(i, j - 1, k) ? j : j - 1][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleU(GetGridCenter(i, j - 1, k - 1) + vector3(stepSize / 2, 0, 0));
			}
			else if (localPos.y >= 0.5 * stepSize && localPos.z < 0.5 * stepSize)
			{
				xValue3 = isVelocityValid('u', i, j, k - 1) ? u[i][j][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleU(GetGridCenter(i, j, k - 1) + vector3(-stepSize / 2, 0, 0));
				xValue4 = isVelocityValid('u', i + 1, j, k - 1) ? u[i + 1][j][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleU(GetGridCenter(i, j, k - 1) + vector3(stepSize / 2, 0, 0));
				yValue1 = isVelocityValid('u', i, j + 1, k) ? u[i][isSolidCell(i, j + 1, k) ? j : j + 1][k] : coarseLayer->SampleU(GetGridCenter(i, j + 1, k) + vector3(-stepSize / 2, 0, 0));
				yValue2 = isVelocityValid('u', i + 1, j + 1, k) ? u[i + 1][isSolidCell(i, j + 1, k) ? j : j + 1][k] : coarseLayer->SampleU(GetGridCenter(i, j + 1, k) + vector3(stepSize / 2, 0, 0));
				yValue3 = isVelocityValid('u', i, j + 1, k - 1) ? u[i][isSolidCell(i, j + 1, k) ? j : j + 1][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleU(GetGridCenter(i, j + 1, k - 1) + vector3(-stepSize / 2, 0, 0));
				yValue4 = isVelocityValid('u', i + 1, j + 1, k - 1) ? u[i + 1][isSolidCell(i, j + 1, k) ? j : j + 1][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleU(GetGridCenter(i, j + 1, k - 1) + vector3(stepSize / 2, 0, 0));
			}

			resultU = TrilinearInterpolate<DOUBLE>(xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, xAlpha, zAlpha, yAlpha);
		}
		if (bPrint)
		{
			(*DecoLogger::GetSingleton()) << "xValue1 :" << xValue1 << "\n";
			(*DecoLogger::GetSingleton()) << "xValue2 :" << xValue2 << "\n";
			(*DecoLogger::GetSingleton()) << "xValue3 :" << xValue3 << "\n";
			(*DecoLogger::GetSingleton()) << "xValue4 :" << xValue4 << "\n";
			(*DecoLogger::GetSingleton()) << "yValue1 :" << yValue1 << "\n";
			(*DecoLogger::GetSingleton()) << "yValue2 :" << yValue2 << "\n";
			(*DecoLogger::GetSingleton()) << "yValue3 :" << yValue3 << "\n";
			(*DecoLogger::GetSingleton()) << "yValue4 :" << yValue4 << "\n";
			(*DecoLogger::GetSingleton()) << "xAlpha :" << xAlpha << "\n";
			(*DecoLogger::GetSingleton()) << "yAlpha :" << yAlpha << "\n";
			(*DecoLogger::GetSingleton()) << "zAlpha :" << zAlpha << "\n";
			(*DecoLogger::GetSingleton()) << "resultU :" << resultU << "\n";
		}
		return resultU;
	}
}
DOUBLE GridField::SampleV(const vector3& pos, GridField* coarseLayer) const
{
	if (this == coarseLayer || !coarseLayer)
		return SampleV(pos);
	else
	{
		DOUBLE resultV;
		vector3 localPos;
		INT numX = numXElement;
		INT numY = numYElement;
		INT numZ = numZElement;
		INT i = 0, j = 0, k = 0;

		localPos = fromWorldToLocalGrid(pos, (UINT&)i, (UINT&)j, (UINT&)k);

		DOUBLE xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, xAlpha, yAlpha, zAlpha; 


		{
			xValue1 = isVelocityValid('v', i, j, k) ? v[i][j][k] : coarseLayer->SampleV(GetGridCenter(i, j, k) + vector3(0, -stepSize / 2, 0));
			xValue2 = isVelocityValid('v', i, j + 1, k) ? v[i][j + 1][k] : coarseLayer->SampleV(GetGridCenter(i, j, k) + vector3(0, stepSize / 2, 0));
			yAlpha = localPos.y / stepSize;
			zAlpha = signedAbs(stepSize / 2 - localPos.z) / stepSize;
			xAlpha = signedAbs(stepSize / 2 - localPos.x) / stepSize;

			if (localPos.x >= 0.5 * stepSize && localPos.z >= 0.5 * stepSize)
			{
				xValue3 = (isVelocityValid('v', i, j, k + 1) ? v[i][j][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleV(GetGridCenter(i, j, k + 1) + vector3(0, -stepSize / 2, 0)));
				xValue4 = (isVelocityValid('v', i, j + 1, k + 1) ? v[i][j + 1][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleV(GetGridCenter(i, j, k + 1) + vector3(0, stepSize / 2, 0)));
				yValue1 = (isVelocityValid('v', i + 1, j, k) ? v[isSolidCell(i + 1, j, k) ? i : i + 1][j][k] : coarseLayer->SampleV(GetGridCenter(i + 1, j, k) + vector3(0, -stepSize / 2, 0)));
				yValue2 = (isVelocityValid('v', i + 1, j + 1, k) ? v[isSolidCell(i + 1, j, k) ? i : i + 1][j + 1][k] : coarseLayer->SampleV(GetGridCenter(i + 1, j, k) + vector3(0, stepSize / 2, 0)));
				yValue3 = (isVelocityValid('v', i + 1, j, k + 1) ? v[isSolidCell(i + 1, j, k) ? i : i + 1][j][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleV(GetGridCenter(i + 1, j, k + 1) + vector3(0, -stepSize / 2,0)));
				yValue4 = (isVelocityValid('v', i + 1, j + 1, k + 1) ? v[isSolidCell(i + 1, j, k) ? i : i + 1][j + 1][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleV(GetGridCenter(i + 1, j , k + 1) + vector3(0, stepSize / 2, 0)));
			}
			else if (localPos.x < 0.5 * stepSize && localPos.z >= 0.5 * stepSize)
			{
				xValue3 = (isVelocityValid('v', i, j, k + 1) ? v[i][j][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleV(GetGridCenter(i, j, k + 1) + vector3(0, -stepSize / 2, 0)));
				xValue4 = (isVelocityValid('v', i, j + 1, k + 1) ? v[i][j + 1][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleV(GetGridCenter(i, j, k + 1) + vector3(0, stepSize / 2, 0)));
				yValue1 = (isVelocityValid('v', i - 1, j, k) ? v[isSolidCell(i - 1, j, k) ? i : i - 1][j][k] : coarseLayer->SampleV(GetGridCenter(i - 1, j, k) + vector3(0, -stepSize / 2, 0)));
				yValue2 = (isVelocityValid('v', i - 1, j + 1, k) ? v[isSolidCell(i - 1, j, k) ? i : i - 1][j + 1][k] : coarseLayer->SampleV(GetGridCenter(i - 1, j, k) + vector3(0, stepSize / 2, 0)));
				yValue3 = (isVelocityValid('v', i - 1, j, k + 1) ? v[isSolidCell(i - 1, j, k) ? i : i - 1][j][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleV(GetGridCenter(i - 1, j, k + 1) + vector3(0, -stepSize / 2, 0)));
				yValue4 = (isVelocityValid('v', i - 1, j + 1, k + 1) ? v[isSolidCell(i - 1, j, k) ? i : i - 1][j + 1][isSolidCell(i, j, k + 1) ? k : k + 1] : coarseLayer->SampleV(GetGridCenter(i - 1, j, k + 1) + vector3(0, stepSize / 2, 0)));
			}
			else if (localPos.x < 0.5 * stepSize && localPos.z < 0.5 * stepSize)
			{
				xValue3 = (isVelocityValid('v', i, j, k - 1) ? v[i][j][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleV(GetGridCenter(i, j, k - 1) + vector3(0, -stepSize / 2, 0)));
				xValue4 = (isVelocityValid('v', i, j + 1, k - 1) ? v[i][j + 1][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleV(GetGridCenter(i, j, k - 1) + vector3(0, stepSize / 2, 0)));
				yValue1 = (isVelocityValid('v', i - 1, j, k) ? v[isSolidCell(i - 1, j, k) ? i : i - 1][j][k] : coarseLayer->SampleV(GetGridCenter(i - 1, j, k) + vector3(0, -stepSize / 2, 0)));
				yValue2 = (isVelocityValid('v', i - 1, j + 1, k) ? v[isSolidCell(i - 1, j, k) ? i : i - 1][j + 1][k] : coarseLayer->SampleV(GetGridCenter(i - 1, j, k) + vector3(0, stepSize / 2, 0)));
				yValue3 = (isVelocityValid('v', i - 1, j, k - 1) ? v[isSolidCell(i - 1, j, k) ? i : i - 1][j][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleV(GetGridCenter(i - 1, j, k - 1) + vector3(0, -stepSize / 2, 0))); 
				yValue4 = (isVelocityValid('v', i - 1, j + 1, k - 1) ? v[isSolidCell(i - 1, j, k) ? i : i - 1][j + 1][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleV(GetGridCenter(i - 1, j, k - 1) + vector3(0, stepSize / 2, 0)));
			}
			else if (localPos.x >= 0.5 * stepSize && localPos.z < 0.5 * stepSize)
			{
				xValue3 = (isVelocityValid('v', i, j, k - 1) ? v[i][j][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleV(GetGridCenter(i, j, k - 1) + vector3(0, -stepSize / 2, 0)));
				xValue4 = (isVelocityValid('v', i, j + 1, k - 1) ? v[i][j + 1][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleV(GetGridCenter(i, j, k - 1) + vector3(0, stepSize / 2, 0)));
				yValue1 = (isVelocityValid('v', i + 1, j, k) ? v[isSolidCell(i + 1, j, k) ? i : i + 1][j][k] : coarseLayer->SampleV(GetGridCenter(i + 1, j, k) + vector3(0, -stepSize / 2, 0)));
				yValue2 = (isVelocityValid('v', i + 1, j + 1, k) ? v[isSolidCell(i + 1, j, k) ? i : i + 1][j + 1][k] : coarseLayer->SampleV(GetGridCenter(i + 1, j, k) + vector3(0, stepSize / 2, 0)));
				yValue3 = (isVelocityValid('v', i + 1, j, k - 1) ? v[isSolidCell(i + 1, j, k) ? i : i + 1][j][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleV(GetGridCenter(i + 1, j, k - 1) + vector3(0, - stepSize / 2, 0)));
				yValue4 = (isVelocityValid('v', i + 1, j + 1, k - 1) ? v[isSolidCell(i + 1, j, k) ? i : i + 1][j + 1][isSolidCell(i, j, k - 1) ? k : k - 1] : coarseLayer->SampleV(GetGridCenter(i + 1, j, k - 1) + vector3(0, stepSize / 2, 0)));
			}
			else 
				assert(0);
			resultV = TrilinearInterpolate<DOUBLE>(xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, yAlpha, zAlpha, xAlpha);
		}
		return resultV;

	}
}
DOUBLE GridField::SampleW(const vector3& pos, GridField* coarseLayer) const
{
	if (this == coarseLayer || !coarseLayer)
		return SampleW(pos);
	else
	{
		DOUBLE resultW;
		vector3 localPos;
		INT numX = numXElement;
		INT numY = numYElement;
		INT numZ = numZElement;
		INT i = 0, j = 0, k = 0;

		localPos = fromWorldToLocalGrid(pos, (UINT&)i, (UINT&)j, (UINT&)k);

		DOUBLE xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, xAlpha, yAlpha, zAlpha; 

		{
			xValue1 = isVelocityValid('w', i, j, k) ? w[i][j][k] : coarseLayer->SampleW(GetGridCenter(i, j, k) + vector3(0, 0, -stepSize / 2));
			xValue2 = isVelocityValid('w', i, j, k + 1) ? w[i][j][k + 1] : coarseLayer->SampleW(GetGridCenter(i, j, k) + vector3(0, 0, stepSize / 2));
			zAlpha = localPos.z / stepSize;
			xAlpha = signedAbs(stepSize / 2 - localPos.x) / stepSize;
			yAlpha = signedAbs(stepSize / 2 - localPos.y) / stepSize;

			if (localPos.y >= 0.5 * stepSize && localPos.x >= 0.5 * stepSize)
			{
				xValue3 = isVelocityValid('w', i + 1, j, k) ? w[isSolidCell(i + 1, j, k) ? i : i + 1][j][k] : coarseLayer->SampleW(GetGridCenter(i + 1,j ,k) + vector3(0, 0, -stepSize / 2));
				xValue4 = isVelocityValid('w', i + 1, j, k + 1) ? w[isSolidCell(i + 1, j, k) ? i : i + 1][j][k + 1] : coarseLayer->SampleW(GetGridCenter(i + 1, j, k) + vector3(0, 0, stepSize / 2));
				yValue1 = isVelocityValid('w', i, j + 1, k) ? w[i][isSolidCell(i, j + 1, k) ? j : j + 1][k] : coarseLayer->SampleW(GetGridCenter(i, j + 1, k) + vector3(0, 0, -stepSize / 2));
				yValue2 = isVelocityValid('w', i, j + 1, k + 1) ? w[i][isSolidCell(i, j + 1, k) ? j : j + 1][k + 1] : coarseLayer->SampleW(GetGridCenter(i, j + 1, k) + vector3(0, 0, stepSize / 2));
				yValue3 = isVelocityValid('w', i + 1, j + 1, k) ? w[isSolidCell(i + 1, j, k) ? i : i + 1][isSolidCell(i, j + 1, k) ? j : j + 1][k] : coarseLayer->SampleW(GetGridCenter(i + 1, j + 1, k) + vector3(0, 0, -stepSize / 2));
				yValue4 = isVelocityValid('w', i + 1, j + 1, k + 1) ? w[isSolidCell(i + 1, j, k) ? i : i + 1][isSolidCell(i, j + 1, k) ? j : j + 1][k + 1] : coarseLayer->SampleW(GetGridCenter(i + 1, j + 1, k) + vector3(0, 0, stepSize / 2));
			}
			else if (localPos.y < 0.5 * stepSize && localPos.x >= 0.5 * stepSize)
			{
				xValue3 = isVelocityValid('w', i + 1, j, k) ? w[isSolidCell(i + 1, j, k) ? i : i + 1][j][k] : coarseLayer->SampleW(GetGridCenter(i + 1,j ,k) + vector3(0, 0, -stepSize / 2));
				xValue4 = isVelocityValid('w', i + 1, j, k + 1) ? w[isSolidCell(i + 1, j, k) ? i : i + 1][j][k + 1] : coarseLayer->SampleW(GetGridCenter(i + 1, j, k) + vector3(0, 0, stepSize / 2));
				yValue1 = isVelocityValid('w', i, j - 1, k) ? w[i][isSolidCell(i, j - 1, k) ? j : j - 1][k] : coarseLayer->SampleW(GetGridCenter(i, j - 1, k) + vector3(0, 0, -stepSize / 2));
				yValue2 = isVelocityValid('w', i, j - 1, k + 1) ? w[i][isSolidCell(i, j - 1, k) ? j : j - 1][k + 1] : coarseLayer->SampleW(GetGridCenter(i, j - 1, k) + vector3(0, 0, stepSize / 2));
				yValue3 = isVelocityValid('w', i + 1, j - 1, k) ? w[isSolidCell(i + 1, j, k) ? i : i + 1][isSolidCell(i, j - 1, k) ? j : j - 1][k] : coarseLayer->SampleW(GetGridCenter(i + 1, j - 1, k) + vector3(0, 0, -stepSize / 2));
				yValue4 = isVelocityValid('w', i + 1, j - 1, k + 1) ? w[isSolidCell(i + 1, j, k) ? i : i + 1][isSolidCell(i, j - 1, k) ? j : j - 1][k + 1] : coarseLayer->SampleW(GetGridCenter(i + 1, j - 1, k) + vector3(0, 0, stepSize / 2));
			}
			else if (localPos.y < 0.5 * stepSize && localPos.x < 0.5 * stepSize)
			{
				xValue3 = isVelocityValid('w', i - 1, j, k) ? w[isSolidCell(i - 1, j, k) ? i : i - 1][j][k] : coarseLayer->SampleW(GetGridCenter(i - 1, j, k) + vector3(0, 0, -stepSize / 2));
				xValue4 = isVelocityValid('w', i - 1, j, k + 1) ? w[isSolidCell(i - 1, j, k) ? i : i - 1][j][k + 1] : coarseLayer->SampleW(GetGridCenter(i - 1, j, k) + vector3(0, 0, stepSize / 2));
				yValue1 = isVelocityValid('w', i, j - 1, k) ? w[i][isSolidCell(i, j - 1, k) ? j : j - 1][k] : coarseLayer->SampleW(GetGridCenter(i, j - 1, k) + vector3(0, 0, -stepSize / 2));
				yValue2 = isVelocityValid('w', i, j - 1, k + 1) ? w[i][isSolidCell(i, j - 1, k) ? j : j - 1][k + 1] : coarseLayer->SampleW(GetGridCenter(i, j - 1, k) + vector3(0, 0, stepSize / 2));
				yValue3 = isVelocityValid('w', i - 1, j - 1, k) ? w[isSolidCell(i - 1, j, k) ? i : i - 1][isSolidCell(i, j - 1, k) ? j : j - 1][k] : coarseLayer->SampleW(GetGridCenter(i - 1, j - 1, k) + vector3(0, 0, -stepSize / 2));
				yValue4 = isVelocityValid('w', i - 1, j - 1, k + 1) ? w[isSolidCell(i - 1, j, k) ? i : i - 1][isSolidCell(i, j - 1, k) ? j : j - 1][k + 1] : coarseLayer->SampleW(GetGridCenter(i - 1, j - 1, k) + vector3(0, 0, stepSize / 2));
			}
			else if (localPos.y >= 0.5 * stepSize && localPos.x < 0.5 * stepSize)
			{
				xValue3 = isVelocityValid('w', i - 1, j, k) ? w[isSolidCell(i - 1, j, k) ? i : i - 1][j][k] : coarseLayer->SampleW(GetGridCenter(i - 1, j, k) + vector3(0, 0, -stepSize / 2));
				xValue4 = isVelocityValid('w', i - 1, j, k + 1) ? w[isSolidCell(i - 1, j, k) ? i : i - 1][j][k + 1] : coarseLayer->SampleW(GetGridCenter(i - 1, j, k) + vector3(0, 0, stepSize / 2));
				yValue1 = isVelocityValid('w', i, j + 1, k) ? w[i][isSolidCell(i, j + 1, k) ? j : j + 1][k] : coarseLayer->SampleW(GetGridCenter(i, j + 1, k) + vector3(0, 0, -stepSize / 2));
				yValue2 = isVelocityValid('w', i, j + 1, k + 1) ? w[i][isSolidCell(i, j + 1, k) ? j : j + 1][k + 1] : coarseLayer->SampleW(GetGridCenter(i, j + 1, k) + vector3(0, 0, stepSize / 2));
				yValue3 = isVelocityValid('w', i - 1, j + 1, k) ? w[isSolidCell(i - 1, j, k) ? i : i - 1][isSolidCell(i, j + 1, k) ? j : j + 1][k] : coarseLayer->SampleW(GetGridCenter(i - 1, j + 1, k) + vector3(0, 0, -stepSize / 2));
				yValue4 = isVelocityValid('w', i - 1, j + 1, k + 1) ? w[isSolidCell(i - 1, j, k) ? i : i - 1][isSolidCell(i, j + 1, k) ? j : j + 1][k + 1] : coarseLayer->SampleW(GetGridCenter(i - 1, j + 1, k) + vector3(0, 0, stepSize / 2));
			}
			else
				assert(0);
			resultW = TrilinearInterpolate<DOUBLE>(xValue1, xValue2, xValue3, xValue4, yValue1, yValue2, yValue3, yValue4, zAlpha, xAlpha, yAlpha);
		}
		return resultW;
	}
}

vector3 GridField::SampleVelocity(const vector3& pos, GridField* coarseLayer) const
{
	vector3 result;
	INT xIndex = 0, yIndex = 0, zIndex = 0;
	const vector3& localPos = fromWorldToLocalGrid(pos, (UINT&)xIndex, (UINT&)yIndex, (UINT&)zIndex);
	result.x = SampleUWithLocalGrid(pos, localPos, xIndex, yIndex, zIndex, coarseLayer);
	result.y = SampleVWithLocalGrid(pos, localPos, xIndex, yIndex, zIndex, coarseLayer);
	result.z = SampleWWithLocalGrid(pos, localPos, xIndex, yIndex, zIndex, coarseLayer);

	return result;
}

vector3 GridField::SampleVelocity(const vector3& pos) const
{
	vector3 result;
	INT xIndex = 0, yIndex = 0, zIndex = 0;
	const vector3& localPos = fromWorldToLocalGrid(pos, (UINT&)xIndex, (UINT&)yIndex, (UINT&)zIndex);
	result.x = SampleUWithLocalGrid(localPos, xIndex, yIndex, zIndex);
	result.y = SampleVWithLocalGrid(localPos, xIndex, yIndex, zIndex);
	result.z = SampleWWithLocalGrid(localPos, xIndex, yIndex, zIndex);

	return result;
}

void GridField::calculateVelocityAtGridCenter()
{
	for (int i = 0; i < numXElement; ++i)
	{
		for (int j = 0; j < numYElement; ++j)
		{
			for (int k = 0; k < numZElement; ++k)
			{
				mVelocity[i][j][k].x = (u[i][j][k] + u[i + 1][j][k]) / 2;
				mVelocity[i][j][k].y = (v[i][j][k] + v[i][j + 1][k]) / 2;
				mVelocity[i][j][k].z = (w[i][j][k] + w[i][j][k + 1]) / 2;
			}
		}
	}
}

double GridField::GetMaxPressureMagnitute() const
{
	double ret = MIN_DOUBLE;
	for (int i = 0; i < numXElement; ++i)
	{
		for (int j = 0; j < numYElement; ++j)
		{
			for (int k = 0; k < numZElement; ++k)
			{
				if (isFluidCell(i, j, k))
				{
					double pressure = p[i][j][k];
					if (pressure > ret)
						ret = pressure;
				}
			}
		}
	}
	return ret;
}
double GridField::GetMinPressureMagnitute() const
{
	double ret = MAX_DOUBLE;
	for (int i = 0; i < numXElement; ++i)
	{
		for (int j = 0; j < numYElement; ++j)
		{
			for (int k = 0; k < numZElement; ++k)
			{
				if (isFluidCell(i, j, k))
				{
					double pressure = p[i][j][k];
					if (pressure < ret)
						ret = pressure;
				}
			}
		}
	}
	return ret;
}


double GridField::GetMaxVorticiyMagnitute() const
{
	double ret = -1;
	for (int i = 0; i < numXElement; ++i)
	{
		for (int j = 0; j < numYElement; ++j)
		{
			for (int k = 0; k < numZElement; ++k)
			{
				double vortMagnitute = mVorticity[i][j][k];
				if (vortMagnitute > ret)
					ret = vortMagnitute;
			}
		}
	}
	return sqrt(ret);
}

bool GridField::MoveGrid(const vector3& offset, ParticleSys* particles, DecoScene* scene)
{
	bool bMoved = false;
	if (abs(offset.x) > abs(scnBox.GetExtent().x) / 2.0)
	{
		moveGridLeftRight(offset.x, particles, scene);
		bMoved = true;
	}
	if (abs(offset.y) > abs(scnBox.GetExtent().y) / 2.0)
	{
		moveGridUpDown(offset.y, particles, scene);
		bMoved = true;
	}
	if (abs(offset.z) > abs(scnBox.GetExtent().z) / 2.0)
	{
		moveGridFrontBack(offset.z, particles, scene);
		bMoved = true;
	}
	if (bMoved)
		*(DecoLogger::GetSingleton()) << "Grid Moved.\n";
	return bMoved;
}

void GridField::moveGridLeftRight(double offset, ParticleSys* particles, DecoScene* scene)
{
	int numGridOffset = static_cast<int>(offset / stepSize);
	if (!numGridOffset) return;
	double realOffset = numGridOffset * stepSize;
	vector3 scnBoxMin = scnBox.GetMin();
	vector3 scnBoxMax = scnBox.GetMax();
	scnBoxMin.x += realOffset;
	scnBoxMax.x += realOffset;

	scene->SetSceneBox(Box(scnBoxMin, scnBoxMax));
	Box expandBox;
	int airMargin = 1;
	DecoConfig::GetSingleton()->GetInt("Solver", "airMargin", airMargin);

	if (numGridOffset > 0)
	{
		for (int i = 0; i < numXElement; ++i)
		{
			for (int j = 0; j < numYElement; ++j)
			{
				for (int k = 0; k < numZElement; ++k)
				{
					if (i + numGridOffset < numXElement)
					{
						u[i][j][k] = u[i + numGridOffset][j][k];
						v[i][j][k] = v[i + numGridOffset][j][k];
						w[i][j][k] = w[i + numGridOffset][j][k];
					}
					else
					{
						u[i][j][k] = 0;
						v[i][j][k] = 0;
						w[i][j][k] = 0;
					}
					u[numXElement][j][k] = 0;
				}
			}
		}
		for (int i = 0; i < numXElement; ++i)
		{
			for (int j = 0; j < numYElement; ++j)
			{
				if (i + numGridOffset < numXElement)
					w[i][j][numZElement] = w[i + numGridOffset][j][numZElement];
				else
					w[i][j][numZElement] = 0;
			}
			for (int k = 0; k < numZElement; ++k)
			{
				if (i + numGridOffset < numXElement)
					v[i][numYElement][k] = v[i + numGridOffset][numYElement][k];
				else
					v[i][numYElement][k] = 0;
			}
		}
		expandBox.Set(vector3(scnBoxMax.x - realOffset - airMargin * stepSize, scnBoxMin.y, scnBoxMin.z), vector3(scnBoxMax.x, scnBoxMax.y, scnBoxMax.z));
	}
	else
	{
		for (int i = numXElement - 1; i >= 0; --i)
		{
			for (int j = 0; j < numYElement; ++j)
			{
				if (i + numGridOffset >= 0)
					w[i][j][numZElement] = w[i + numGridOffset][j][numZElement];
				else
					w[i][j][numZElement] = 0;
			}
			for (int k = 0; k < numZElement; ++k)
			{
				if (i + numGridOffset >= 0)
					v[i][numYElement][k] = v[i + numGridOffset][numYElement][k];
				else
					v[i][numYElement][k] = 0;
			}
		}

		for (int i = numXElement - 1; i >= 0; --i)
		{
			for (int j = 0; j < numYElement; ++j)
			{
				for (int k = 0; k < numZElement; ++k)
				{
					if (i == numXElement - 1)
						u[numXElement][j][k] = u[numXElement + numGridOffset][j][k];

					if (i + numGridOffset >= 0)
					{
						u[i][j][k] = u[i + numGridOffset][j][k];
						v[i][j][k] = v[i + numGridOffset][j][k];
						w[i][j][k] = w[i + numGridOffset][j][k];
					}
					else
					{
						u[i][j][k] = 0;
						v[i][j][k] = 0;
						w[i][j][k] = 0;
					}
				}
			}
		}
		expandBox.Set(vector3(scnBoxMin.x, scnBoxMin.y, scnBoxMin.z), vector3(scnBoxMin.x - realOffset + airMargin * stepSize, scnBoxMax.y, scnBoxMax.z));
	}
	particles->EraseParticlesOutsideBox(scnBox);
	scene->VoxelizeSolid();
	scene->solver.UpdateFluidGridType(false);
	particles->FillParticlesInBox(expandBox);
}

void GridField::moveGridUpDown(double offset, ParticleSys* particles, DecoScene* scene)
{
	int numGridOffset = static_cast<int>(offset / stepSize);
	if (!numGridOffset) return;
	double realOffset = numGridOffset * stepSize;
	vector3 scnBoxMin = scnBox.GetMin();
	vector3 scnBoxMax = scnBox.GetMax();
	scnBoxMin.y += realOffset;
	scnBoxMax.y += realOffset;

	scene->SetSceneBox(Box(scnBoxMin, scnBoxMax));

	Box expandBox;
	int airMargin = 1;
	DecoConfig::GetSingleton()->GetInt("Solver", "airMargin", airMargin);

	if (numGridOffset > 0)
	{
		expandBox.Set(vector3(scnBoxMin.x, scnBoxMax.y - realOffset - airMargin * stepSize, scnBoxMin.z), vector3(scnBoxMax.x, scnBoxMax.y, scnBoxMax.z));
		for (int i = 0; i < numXElement; ++i)
		{
			for (int j = 0; j < numYElement; ++j)
			{
				for (int k = 0; k < numZElement; ++k)
				{
					if (j + numGridOffset < numYElement)
					{
						u[i][j][k] = u[i][j + numGridOffset][k];
						v[i][j][k] = v[i][j + numGridOffset][k];
						w[i][j][k] = w[i][j + numGridOffset][k];
					}
					else
					{
						u[i][j][k] = 0;
						v[i][j][k] = 0;
						w[i][j][k] = 0;
					}
					v[i][numYElement][k] = 0;
				}
			}
		}
		for (int j = 0; j < numYElement; ++j)
		{
			for (int i = 0; i < numXElement; ++i)
			{
				if (j + numGridOffset < numYElement)
					w[i][j][numZElement] = w[i][j + numGridOffset][numZElement];
				else
					w[i][j][numZElement] = 0;
			}
			for (int k = 0; k < numZElement; ++k)
			{
				if (j + numGridOffset < numYElement)
					u[numXElement][j][k] = u[numXElement][j + numGridOffset][k];
				else
					u[numXElement][j][k] = 0;
			}
		}
	}
	else
	{
		expandBox.Set(vector3(scnBoxMin.x, scnBoxMin.y, scnBoxMin.z), vector3(scnBoxMax.x, scnBoxMin.y - realOffset + airMargin * stepSize, scnBoxMax.z));
		for (int j = numYElement - 1; j >= 0; --j)
		{
			for (int i = 0; i < numXElement; ++i)
			{
				if (j + numGridOffset >= 0)
					w[i][j][numZElement] = w[i][j + numGridOffset][numZElement];
				else
					w[i][j][numZElement] = 0;
			}
			for (int k = 0; k < numZElement; ++k)
			{
				if (j + numGridOffset >= 0)
					u[numXElement][j][k] = u[numXElement][j + numGridOffset][k];
				else
					u[numXElement][j][k] = 0;
			}
		}

		for (int i = 0; i < numXElement; ++i)
		{
			for (int j = numYElement - 1; j >= 0; --j)
			{
				for (int k = 0; k < numZElement; ++k)
				{
					if (j == numYElement - 1)
						v[i][numYElement][k] = v[i][numYElement + numGridOffset][k];

					if (j + numGridOffset >= 0)
					{
						u[i][j][k] = u[i][j + numGridOffset][k];
						v[i][j][k] = v[i][j + numGridOffset][k];
						w[i][j][k] = w[i][j + numGridOffset][k];
					}
					else
					{
						u[i][j][k] = 0;
						v[i][j][k] = 0;
						w[i][j][k] = 0;
					}
				}
			}
		}
	}
	scene->VoxelizeSolid();
	scene->solver.UpdateFluidGridType(false);
	particles->EraseParticlesOutsideBox(scnBox);
	particles->FillParticlesInBox(expandBox);
}

void GridField::moveGridFrontBack(double offset, ParticleSys* particles, DecoScene* scene)
{
	int numGridOffset = static_cast<int>(offset / stepSize);
	if (!numGridOffset) return;
	double realOffset = numGridOffset * stepSize;
	vector3 scnBoxMin = scnBox.GetMin();
	vector3 scnBoxMax = scnBox.GetMax();
	scnBoxMin.z += realOffset;
	scnBoxMax.z += realOffset;

	scene->SetSceneBox(Box(scnBoxMin, scnBoxMax));
	int airMargin = 1;
	DecoConfig::GetSingleton()->GetInt("Solver", "airMargin", airMargin);

	Box expandBox;

	if (numGridOffset > 0)
	{
		for (int i = 0; i < numXElement; ++i)
		{
			for (int j = 0; j < numYElement; ++j)
			{
				for (int k = 0; k < numZElement; ++k)
				{
					if (k + numGridOffset < numZElement)
					{
						u[i][j][k] = u[i][j][k + numGridOffset];
						v[i][j][k] = v[i][j][k + numGridOffset];
						w[i][j][k] = w[i][j][k + numGridOffset];
					}
					else
					{
						u[i][j][k] = 0;
						v[i][j][k] = 0;
						w[i][j][k] = 0;
					}
					w[i][j][numZElement] = 0;
				}
			}
		}
		for (int k = 0; k < numZElement; ++k)
		{
			for (int j = 0; j < numYElement; ++j)
			{
				if (k + numGridOffset < numZElement)
					u[numXElement][j][k] = u[numXElement][j][k + numGridOffset];
				else
					u[numXElement][j][k] = 0;
			}
			for (int i = 0; i < numXElement; ++i)
			{
				if (k + numGridOffset < numZElement)
					v[i][numYElement][k] = v[i][numYElement][k + numGridOffset];
				else
					v[i][numYElement][k] = 0;
			}
		}
		expandBox.Set(vector3(scnBoxMin.x, scnBoxMin.y, scnBoxMax.z - realOffset - airMargin * stepSize), vector3(scnBoxMax.x, scnBoxMax.y, scnBoxMax.z));
	}
	else
	{
		for (int k = numZElement - 1; k >= 0 ; --k)
		{
			for (int j = 0; j < numYElement; ++j)
			{
				if (k + numGridOffset >= 0)
					u[numXElement][j][k] = u[numXElement][j][k + numGridOffset];
				else
					u[numXElement][j][k] = 0;
			}
			for (int i = 0; i < numXElement; ++i)
			{
				if (k + numGridOffset >= 0)
					v[i][numYElement][k] = v[i][numYElement][k + numGridOffset];
				else
					v[i][numYElement][k] = 0;
			}
		}

		for (int i = 0; i < numXElement; ++i)
		{
			for (int j = 0; j < numYElement; ++j)
			{
				for (int k = numZElement - 1; k >= 0 ; --k)
				{
					if (k == numZElement - 1)
						w[i][j][numZElement] = u[i][j][numZElement + numGridOffset];

					if (k + numGridOffset >= 0)
					{
						u[i][j][k] = u[i][j][k + numGridOffset];
						v[i][j][k] = v[i][j][k + numGridOffset];
						w[i][j][k] = w[i][j][k + numGridOffset];
					}
					else
					{
						u[i][j][k] = 0;
						v[i][j][k] = 0;
						w[i][j][k] = 0;
					}
				}
			}
		}
		expandBox.Set(vector3(scnBoxMin.x, scnBoxMin.y, scnBoxMin.z), vector3(scnBoxMax.x, scnBoxMax.y, scnBoxMin.z - realOffset + airMargin * stepSize));
	}
	scene->VoxelizeSolid();
	scene->solver.UpdateFluidGridType(false);
	particles->EraseParticlesOutsideBox(scnBox);
	particles->FillParticlesInBox(expandBox);
}


void GridField::CalculateVorticity()
{
	calculateVelocityAtGridCenter();
	double dwdy = 0;
	double dvdz = 0;
	double dudz = 0;
	double dwdx = 0;
	double dvdx = 0;
	double dudy = 0;
	for (int i = 0; i < numXElement; ++i)
	{
		for (int j = 0; j < numYElement; ++j)
		{
			for (int k = 0; k < numZElement; ++k)
			{
				if (i + 1 >= numXElement)
				{
					dwdx = (w[i][j][k] - w[i - 1][j][k]) / (stepSize);
					dvdx = (v[i][j][k] - v[i - 1][j][k]) / (stepSize);

				}
				else if (i - 1 < 0)
				{
					dwdx = (w[i + 1][j][k] - w[i][j][k]) / (stepSize);
					dvdx = (v[i + 1][j][k] - v[i][j][k]) / (stepSize);

				}
				else
				{
					dwdx = (w[i + 1][j][k] - w[i - 1][j][k]) / (2 * stepSize);
					dvdx = (v[i + 1][j][k] - v[i - 1][j][k]) / (2 * stepSize);

				}
				if (j + 1 >= numYElement)
				{
					dwdy = (w[i][j][k] - w[i][j - 1][k]) / (stepSize);
					dudy = (u[i][j][k] - u[i][j - 1][k]) / (stepSize);

				}
				else if (j - 1 < 0)
				{
					dwdy = (w[i][j + 1][k] - w[i][j][k]) / (stepSize);
					dudy = (u[i][j + 1][k] - u[i][j][k]) / (stepSize);

				}
				else
				{
					dwdy = (w[i][j + 1][k] - w[i][j - 1][k]) / (2 * stepSize);
					dudy = (u[i][j + 1][k] - u[i][j - 1][k]) / (2 * stepSize);

				}
				if (k + 1 >= numZElement)
				{
					dvdz = (v[i][j][k] - v[i][j][k - 1]) / (stepSize);
					dudz = (u[i][j][k] - u[i][j][k - 1]) / (stepSize);

				}
				else if (k - 1 < 0)
				{
					dvdz = (v[i][j][k + 1] - v[i][j][k]) / (stepSize);
					dudz = (u[i][j][k + 1] - u[i][j][k]) / (stepSize);

				}
				else
				{
					dvdz = (v[i][j][k + 1] - v[i][j][k - 1]) / (2 * stepSize);
					dudz = (u[i][j][k + 1] - u[i][j][k - 1]) / (2 * stepSize);

				}
				vector3 vorticity(dwdy - dvdz, dudz - dwdx, dvdx - dudy);
				mVorticity[i][j][k] = vorticity.length();
			}
		}
	}

}


void GridField::DumpSolid()
{
	FILE* fp = fopen("SolidGridFile.txt", "w");
	UINT numElement = numXElement * numYElement * numZElement;
	if (!numZElement) return;
	for (INT xIndex = 0; xIndex < numXElement; xIndex++)
	{
		for (INT yIndex = 0; yIndex < numYElement; yIndex++)
		{
			for (INT zIndex = 0; zIndex < numZElement; zIndex++)
			{
				if (type[xIndex][yIndex][zIndex] != GT_Solid)
					continue;
				vector3 center = GetGridCenter(xIndex, yIndex, zIndex);
				double sideLength = 2 * stepSize;
				fprintf(fp, "(%f,%f,%f), %f\n", center.x, center.y, center.z, sideLength);
			}
		}
	}
	fclose(fp);

}


DecoArchive& operator<< (DecoArchive& Ar, const GridField& field)
{
	Ar << field.scnBox;
	Ar << field.numXElement;
	Ar << field.numYElement;
	Ar << field.numZElement;
	Ar << field.stepSize;
	
	INT numX = field.numXElement;
	INT numY = field.numYElement;
	INT numZ = field.numZElement;
	for (INT i = 0; i < numX + 1; i++)
	{
		for (INT j = 0; j < numY; j++)
		{
			for (INT k = 0; k < numZ; k++)
			{
				Ar << field.u[i][j][k];
			}
		}
	}
	for (INT i = 0; i < numX; i++)
	{
		for (INT j = 0; j < numY + 1; j++)
		{
			for (INT k = 0; k < numZ; k++)
			{
				Ar << field.v[i][j][k];
			}
		}
	}
	for (INT i = 0; i < numX; i++)
	{
		for (INT j = 0; j < numY; j++)
		{
			for (INT k = 0; k < numZ + 1; k++)
			{
				Ar << field.w[i][j][k];

			}
		}
	}

	for (INT i = 0; i < numX; i++)
	{
		for (INT j = 0; j < numY; j++)
		{
			for (INT k = 0; k < numZ; k++)
			{
				Ar << field.inDomain[i][j][k];
				Ar << field.prevInDomain[i][j][k];
				Ar << field.type[i][j][k];
				//Ar << field.p[i][j][k];
				Ar << field.phi[i][j][k];
				//Ar << field.tmpInDomain[i][j][k];
				//Ar << field.pTheta[i][j][k];
				//Ar << field.mTheta[i][j][k];
			}
		}
	}
	Ar.Flush();

	return Ar;
}
DecoArchive& operator>> (DecoArchive& Ar, GridField& field)
{	
	Ar >> field.scnBox;
	Ar >> field.numXElement;
	Ar >> field.numYElement;
	Ar >> field.numZElement;
	Ar >> field.stepSize;


	INT numX = field.numXElement;
	INT numY = field.numYElement;
	INT numZ = field.numZElement;
	field.SetNumElement(numX, numY, numZ);

	for (INT i = 0; i < numX + 1; i++)
	{
		for (INT j = 0; j < numY; j++)
		{
			for (INT k = 0; k < numZ; k++)
			{
				Ar >> field.u[i][j][k];
			}
		}
	}
	for (INT i = 0; i < numX; i++)
	{
		for (INT j = 0; j < numY + 1; j++)
		{
			for (INT k = 0; k < numZ; k++)
			{
				Ar >> field.v[i][j][k];
			}
		}
	}
	for (INT i = 0; i < numX; i++)
	{
		for (INT j = 0; j < numY; j++)
		{
			for (INT k = 0; k < numZ + 1; k++)
			{
				Ar >> field.w[i][j][k];
			}

		}
	}
	for (INT i = 0; i < numX; i++)
	{
		for (INT j = 0; j < numY; j++)
		{
			for (INT k = 0; k < numZ; k++)
			{
				Ar >> field.inDomain[i][j][k];
				Ar >> field.prevInDomain[i][j][k];
				Ar >> (INT&)field.type[i][j][k];

				//Ar >> field.p[i][j][k];
				Ar >> field.phi[i][j][k];
				//Ar >> field.tmpInDomain[i][j][k];
				//Ar >> field.pTheta[i][j][k];
				//Ar >> field.mTheta[i][j][k];
			}
		}
	}
	return Ar;
}

//  [4/23/2008 HATEVOL] ADDED start
inline
INT GridField::normalizeIndex(INT tmpIndex, INT numElement) const
{
	//cast tmpIndex ranging from 0 to numElement-1 to return value 0
	//negative tmpIndex to negative return value while tmpIndex larger than or equal to numElement to positive return value
	//without branch
	INT numElement5 = numElement * 5;
	return ( (tmpIndex * 10 - (numElement5 - 5)) / (numElement5) );
}
//  [4/23/2008 HATEVOL] ADDED end
